;\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
FUNCTION AREA_AV,data

;Check to see the number of DIMENSIONS of data

;LOOP over the first two dimensions 
FOR I=0,1 DO BEGIN
	sav=size(data)
	data=TOTAL(TEMPORARY(data),1)/sav(1)
ENDFOR

av=data
RETURN,av
END
;\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
PRO CLIMATE,time,dt,ac,nh,clim
	s=size(ac)
	pi=4d*ATAN(1d)
	IF (s(0) gt 2) THEN BEGIN
		nx=s(1) & ny=s(2)
		IF (s(0) gt 3) THEN nz=s(3) ELSE nz=1
	ENDIF ELSE BEGIN
		IF s(0) gt 1 THEN BEGIN
		ny=s(1) & nx=1 & nz=1
		ac=REFORM(ac,nx,ny,nz,s(2))
		ENDIF ELSE BEGIN
		nx=1 & ny=1 & nz=1
		; no need to reform since coming in 1,1,nh+1
		ENDELSE 
	ENDELSE
	nt=N_ELEMENTS(time)
	nt=LONG(nt)
	
	wt=2D*pi*time/(dt*365.25D)
	;HLP,nx,ny,nh,ac
	annual=REFORM(ac,nx*ny*nz,nh+1)
	; add the cmean
	clim = REBIN(annual(*,0),nx*ny*nz,nt)

	FOR j=1,nh,2 DO BEGIN
		fact=DOUBLE((j+1)/2)
		clim=TEMPORARY(clim)+COS(fact*wt)##annual(*,j)
		clim=TEMPORARY(clim)+SIN(fact*wt)##annual(*,j+1)
	ENDFOR
	
	

	clim=REFORM(TEMPORARY(clim),nx,ny,nz,nt)
	clim=REFORM(clim)
	RETURN
End
;//////////////////////////////////////////////////////////
PRO CLIMSER,t,dt,annual_coef,nh,clim
	; t contains the time values
	; dt contains the time units of 1 day
	; ex: if t is in hours,dt=24.0D0
	; annual_coef are the harmonic coefficients
	;	at grid points
	; nh is the number of harmonics*2
	; clim is the output harmonic reconstruction
	time=t
	s=size(annual_coef)
	pi=4d*ATAN(1d)
	nt=N_ELEMENTS(time)
	IF nt eq 1 THEN time=[time]
	; need to determine if annual_coef is 4d,3d,2d, or 1d
	
	CASE s(0) OF
		4:BEGIN
			nx=s(1) & ny=s(2) & nz=s(3)
			clim=FLTARR(nx,ny,nz,nt,/noz)
			
			FOR j=0,ny-1 DO BEGIN
			FOR k=0,nz-1 DO BEGIN
			FOR i=0,nx-1 DO BEGIN
			 CLIMATE,time,dt,annual_coef(i,j,k,*),nh,hclim
			 clim(i,j,k,*)=hclim
			ENDFOR;k
			ENDFOR;j
			ENDFOR;i
		END
	        3:BEGIN
			nx=s(1) & ny=s(2)
			clim=FLTARR(nx,ny,nt,/noz)
			FOR j=0,ny-1 DO BEGIN
			FOR i=0,nx-1 DO BEGIN
		 	 CLIMATE,time,dt,annual_coef(i,j,*),nh,hclim
		 	 clim(i,j,*)=hclim
			ENDFOR;i
			ENDFOR;j
		END		
		ELSE:BEGIN
			IF s(0) gt 1 THEN BEGIN
				ny=s(1) & nx=1
				annual_coef=REFORM(annual_coef,nx,ny,s(2))
				clim=FLTARR(nx,ny,nt,/noz)
				FOR i=0,nx-1 DO BEGIN
				FOR j=0,ny-1 DO BEGIN
					CLIMATE,time,dt,annual_coef(i,j,*),nh,hclim
					clim(i,j,*)=hclim
				ENDFOR
				ENDFOR
			ENDIF ELSE BEGIN
				nx=1 & ny=1
				annual_coef=REFORM(annual_coef,1,1,s(1))
				CLIMATE,time,dt,annual_coef,nh,hclim
				clim=hclim

			ENDELSE 
		END

	ENDCASE

	
RETURN
END
;+
; NAME:
;    HARMFIT
;
; PURPOSE:
;    Harmonic series fit to a function.
;
; CATEGORY:
;    E2 - Curve and Surface Fitting.
;
; CALLING SEQUENCE:
;   YFIT = HARMFIT(X,Y,PERIOD,HARMONIC,ZERO=zero)
;
; INPUTS:
;
;    X:  A row vector of independent variables.
;
;    Y:  A row vector of dependent variable, the same length as x.
;
;    PERIOD: The period of the 1st harmonic cos & sin function
;
;    HARMONIC:  The harmonic to fit up to. i.e. HARMONIC=1 fits
;      up to the 1st harmonic, HARMONIC=2 fits up to the 2nd harmonic
;      (nth harmonic = PERIOD/n)
;
;
; OUTPUTS:
;    YFIT: the vector of the fit.
;
;
; OPTIONAL KEYWORD PARAMETERS:
;
;    ZERO:  A flag to remove the mean (0th harmonic) of the function,
;    or leave the mean in (zero=0). Default is zero=1. If zero is
;    declared beforehand with zero=1 then a call with 'zero=zero'
;    will return the coefficients of the fit in variable zero
;    in the order:
;       [1st harm cos, 1st sin, 2nd cos, 2nd sin, ..., mean]
; 
;    COEFF: if set, then Result is just the coefficients instead of
;           the actual fit. This saves one # operation.
;           Actual fit can be reconstructed by 
;               M + C1*Cos(x) + S1*Sin(x) + ...
;           where M is mean, and C1,S1,C2,S2... are the cos & sin coeffs.
;
;    XFIT: this is an optional X array to use for the returned YFIT.
;          YFIT will then be calculated at these X values rather
;          than the input X values. Useful if you only want to fit
;          using a certain range, but want to remove this fit
;          from a different range.
;
;    DOUBLE: use double precision for calculations.
;
;
; RESTRICTIONS:
; Uses NR_LUDCMP and NR_LUBKSB.
; x and y must have same # of elements.
;
; MODIFICATION HISTORY:
; Written, E. Leuliette & C. Torrence, April 15, 1993.
; Modified & Fixed 4/4/94 to use NR_SVD & NR_SVBKSB (C.T.)
;
;-

FUNCTION HARMFIT,X,Y,PERIOD,HARMONIC, $
	ZERO=zero,COEFF=coeff,XFIT=xfit,DOUBLE=double

	ON_ERROR,2  ;RETURN TO CALLER IF ERROR
	N = N_ELEMENTS(x) 	;SIZE
	IF N NE N_ELEMENTS(y) THEN Begin  ;SAME # OF DATA POINTS.
		message, 'X and Y must have same # of elements.'
		Return,undefined
	ENDIF
	double = KEYWORD_SET(double)
	IF N_ELEMENTS(zero) LE 0 THEN zero = 1
	IF N_ELEMENTS(zero) GT 1 THEN zero = 1

	Pi = 3.1415926D
	omegatime = 2D*Pi/period*DOUBLE(REFORM(x))
	IF NOT double THEN omegatime=FLOAT(omegatime)
	terms = (harmonic*2) + (zero EQ 1)            ;# of terms in fit
	comp = FLTARR(terms,N) + 1.
	IF double THEN comp = DOUBLE(comp)

	FOR i=0,harmonic-1 DO BEGIN
		comp(i*2,*) = COS(omegatime*(i+1))
		comp(i*2+1,*) = SIN(omegatime*(i+1))
	ENDFOR
	good = WHERE(TOTAL(comp^2,2) GT 1./terms,terms)
	comp = comp(good,*)

	matrix = comp # TRANSPOSE(comp)
	B = comp # REFORM(y)
	NR_SVD,matrix,indx,U,V,DOUBLE=double
	B = NR_SVBKSB(U,indx,V,B,DOUBLE=double)

	zero = B  ;return coefficients
	IF (N_ELEMENTS(xfit) GT 0) THEN BEGIN
		omegatime = 2D*Pi/period*DOUBLE(REFORM(xfit))
		IF NOT double THEN omegatime=FLOAT(omegatime)
		comp = FLTARR(terms,N_ELEMENTS(xfit)) + 1.
		FOR i=0,harmonic-1 DO BEGIN
			comp(i*2,*) = COS(omegatime*(i+1))
			comp(i*2+1,*) = SIN(omegatime*(i+1))
		ENDFOR
	ENDIF
	IF KEYWORD_SET(coeff) THEN RETURN,zero ELSE RETURN,REFORM(B # comp)
END
;-------------------------------------------------------------
;+
; NAME:
;       HLP
;
; PURPOSE:
;       Variant of HELP.  Gives array min, max.
;
; CALLING SEQUENCE:
;       HLP, a1, [a2, ..., a9]
;
; INPUTS:
;       a1, [...] = input variables.    in
;
; KEYWORD PARAMETERS:
;       mean = print out mean
;       miss = missing value to exclude from min,max,& mean
;       names = the variable names (string array)
;
; OUTPUTS:
;   gives variable type 
;   for arrays, gives max & min with indices in parentheses
;
; Example:
;
;   IDL> x = [1,5,10,6,4,0]
;   IDL> hlp,x
;      1   INTARR(6) Min(5) = 0, Max(2) = 10
;   IDL> hlp,x,/m,NAMES='X'
;      X   INTARR(6) Min(5)= 0, Max(2)= 10 Mean= 4.33333
;
; MODIFICATION HISTORY:
;       Written R. Sterner, 11 Dec, 1989
;       Modified C. Torrence 14 Apr, 1994
;       Modified C. Torrence 14 Nov, 1995 to use INDICES.pro
;       Modified C. Torrence June 1996 to give Structure information
;-
;-------------------------------------------------------------
 
PRO HLP, a1,a2,a3,a4,a5,a6,a7,a8,a9,MEANN=meann,MISS=miss,NAMES=names

ON_ERROR,2
np = N_PARAMS(0); How many args?
IF (np LT 1) THEN BEGIN
	MESSAGE,/INFO,' Variant of HELP.  Gives array min, max.'
	PRINT,' hlp, a1, [a2, ..., a9]     where a1, [...] = input variables.'
	RETURN
ENDIF

meana=KEYWORD_SET(meann)
missing=N_ELEMENTS(miss) EQ 1
IF NOT KEYWORD_SET(names) THEN names=ROUND_ANY(INDGEN(np)+1)

typ=['UND','BYT','INT','LON','FLT','DBL','CMPLX','STR','STRUC','DCPLX','PTR']

FOR k = 1, np DO BEGIN; Loop through all args.
	kstr=STRTRIM(k,2)
	tmp = EXECUTE('i = N_ELEMENTS(a'+kstr+')') ; Test if ai defined.
	IF i EQ 0 THEN BEGIN
		PRINT,'   '+names(k-1),'   Undefined'
		GOTO, undef
	ENDIF

	tmp = EXECUTE('a = TEMPORARY(a'+kstr+')') ; Put ai in a.
 
	sz = SIZE(a)
	f = FIX(sz(0))		;ISARRAY(a)
	t = sz(f+1)			;DATATYPE(a,3) get variable type

	IF (f GE 1) THEN BEGIN
		aa = 'ARR('
		FOR i = 1, f DO BEGIN
			aa = aa + STRTRIM(sz(i),2)
			IF i LT f THEN aa = aa + ','
		ENDFOR
		aa = aa+')'
		IF (t LT 6) THEN BEGIN
			mni=1
			mxi=1
			IF NOT missing THEN BEGIN
				mn = MIN(a,mni)*1
				mx = MAX(a,mxi)*1
			ENDIF ELSE BEGIN
				mn=WHERE(a NE miss,ngood)
				IF ngood GT 0 THEN BEGIN
					mn = MIN(a(mn),max=mx)*1
					mx=mx*1
					mni=(WHERE(a EQ mn))(0)
					mxi=(WHERE(a EQ mx))(0)
				ENDIF ELSE BEGIN
					mn = miss & mx=miss & mni=0 & mxi=0
				ENDELSE
			ENDELSE
			mni=INDICES(sz,mni)
			mxi=INDICES(sz,mxi)
			mni=STRING(mni,FORMAT='('+STRTRIM(sz(0))+'(I,:,","))')
			mxi=STRING(mxi,FORMAT='('+STRTRIM(sz(0))+'(I,:,","))')
			mni= '('+STRCOMPRESS(mni,/REM)+')'
			mxi= '('+STRCOMPRESS(mxi,/REM)+')'
			mn = STRTRIM(mn,2)
			mx = STRTRIM(mx,2)
			aa=aa+', MIN'+mni+'= '+mn+', MAX'+mxi+'= '+mx
			n = sz(f+2)
			IF meana THEN BEGIN
				aa = aa + ' MEAN= '
				IF NOT missing THEN aa = aa +STRTRIM(TOTAL(a)/n,2) $
				ELSE aa = aa +STRTRIM(TOTAL(a(WHERE(a NE miss,n)))/n,2)
			ENDIF
		ENDIF ELSE BEGIN
			CASE (t) OF
				8: BEGIN
					HELP,a
					HELP,a,/STRUCTURE
					GOTO,skip
				   END
				6: BEGIN
					PRINT,names(k-1),typ(t),aa,FORMAT='(A6,1X,A5,A)'
					HLP,FLOAT(a),IMAGINARY(a),NAMES=['real=','imag=']
					GOTO,skip
				   END
				9: BEGIN
					PRINT,names(k-1),typ(t),aa,FORMAT='(A6,1X,A5,A)'
					HLP,DOUBLE(a),IMAGINARY(a),NAMES=['real=','imag=']
					GOTO,skip
				   END
				ELSE:
			ENDCASE
		ENDELSE
	ENDIF ELSE BEGIN
		IF (typ(t) EQ 'PTR') THEN BEGIN
			PRINT,names(k-1),typ(t),'-->',FORMAT='(A6,1X,A5,A,$)'
			HLP,*a,MEANN=meann,MISS=miss,NAMES=' '
			GOTO,skip
		ENDIF ELSE aa = ' = '+STRTRIM(a,2)
	ENDELSE
	PRINT,names(k-1),typ(t),aa,FORMAT='(A6,1X,A5,A)'
skip:
	tmp = EXECUTE('a'+kstr+' = TEMPORARY(a)') ; Put a back in ai.
undef:
ENDFOR  ; k
 
RETURN
END
;-------------------------------------------------------------
;+
; NAME:
;       INDICES
;
; PURPOSE:
;       Given the SIZE(array) and 1D index numbers, returns
;       the actual indices (x,y,z,...)
;
; CALLING SEQUENCE:
;       ind=INDICES(SIZE(A),N)
;
; INPUTS:
;       SIZE(A) = the vector returned from the SIZE function
;                 (this avoids passing in the entire array)
;       N = the array index (or indices) to convert to true indices
;           (this probably came from a WHERE command)
;
; OUTPUTS:
;       ind = the actual indices of that index.
;       If N is an array (say from WHERE) then the output will be
;       a 2D array of size (ndim,nN) where ndim is the # of dims of A,
;       and nN is the # of N's
;
; EXAMPLE:
;   
;   IDL> A=FLTARR(144,73,14)
;   IDL> A(30,17,3)= 25.3
;   IDL> A(45,12,9)= 35.2
;   IDL> w=WHERE(A GT 0.0)
;   IDL> PRINT,w
;          34014       96381      <--- this tells you nothing!
;   IDL> PRINT,INDICES(SIZE(A),w)
;             30          17           3   <--- this tells you something!
;             45          12           9
;
; BONUS FEATURE:
;       This also does conversion of a # to any base!
;       CALLING SEQUENCE:    ind=INDICES(base,N)
;       OUTPUT: the representation of N in that base.
;
;       EXAMPLE:
;         IDL> PRINT,INDICES(2,12)
;                 1           1           0           0
;
; MODIFICATION HISTORY:
;       Written C.Torrence & G.Compo Nov 1995
;-
;-------------------------------------------------------------


FUNCTION indices,siz1,i1,D=d

	IF (N_ELEMENTS(siz1) GT 32) THEN MESSAGE, $
		'Too many dimensions. Input SIZE(array), not array itself'
	base=-1
	IF (N_ELEMENTS(siz1) EQ 1) THEN BEGIN   ; base conversion
		base=siz1
		siz1=LONG((ALOG(MAX(i1))/ALOG(siz1))+1.0)   ; max # of digits
		IF (N_ELEMENTS(d) GT 0) THEN siz1=d
		siz1=[siz1,LONARR(siz1)+base,3,LONG(base)^siz1]
	ENDIF
	ndim=siz1(0)
	nmax=siz1(ndim+2) 
	ni=N_ELEMENTS(i1)
	
	bad=WHERE(i1 GT nmax,nbad)
	IF nbad GT 0 THEN MESSAGE,'Indices are outside array bounds (N='+ $
			STRCOMPRESS(nmax)+')'
	IF (ndim LE 1) THEN RETURN,i1    ; just a simple vector array
		
	siz=siz1(1:ndim)
	IF base EQ 2 THEN index1=BYTARR(ndim,ni) ELSE index1=LONARR(ndim,ni)
	mult=LONARR(ndim)+siz(0)
	FOR d=1,ndim-2 DO mult(d)=mult(d-1)*siz(d)
	
	FOR i=0,ni-1 DO BEGIN
		mult(ndim-1)=i1(i)+1
		temp=i1(i)
		FOR d=ndim-1,1,-1 DO BEGIN
		   temp=temp MOD mult(d)
		   index1(d,i)=temp/mult(d-1)
		ENDFOR
		index1(0,i)=temp MOD mult(0)
	ENDFOR
	IF base NE -1 THEN index1=REVERSE(index1)
	RETURN,index1
	
END
;//////////////////////////////////////////////////////////////////////////////
FUNCTION REMIX,A,B,NOMATCH=nomatch
;
; NAME:
;	REMIX
;
; PURPOSE:
;	Finds the ordering of one array in terms of another.
;
;	Essentially, this is a generalization of WHERE for two arrays
;	  rather than an array and a scalar.
;
; CATEGORY:
;
; CALLING SEQUENCE:
;	W = REMIX(A,B,[NOMATCH])
;
; INPUTS:
;	A: The first array.
;
;	B: The second array.
;
; OUTPUTS:
;	W: An array of index values such that A=B(W).
;
; OPTIONAL OUTPUT PARAMETERS:
;	NOMATCH: Values in A which do not match any in B.
;
; MODIFICATION HISTORY:
;	Written by: C. Torrence  and E. Leuliette.
;
;

ON_ERROR,2
na=N_ELEMENTS(A)
nb=N_ELEMENTS(B)
nuniq=N_ELEMENTS(UNIQ(B,SORT(B)))

IF (nuniq NE nb) THEN $
	MESSAGE,'b has '+STRCOMPRESS(nb-nuniq,/rem)+' non-unique elements'
IF (na LT 1) THEN MESSAGE,'A has no elements'

c=[-1]

FOR s=0L,na-1 DO c=[c,WHERE(A(s) EQ B)]

c=c(1:*)

nomatch=WHERE(c EQ -1,nno)
IF nno GT 0 THEN nomatch=nomatch

IF (nno LT na) THEN RETURN,c(WHERE(c GT -1)) ELSE RETURN,-1

END

;///////////////////////////////////////////////////////////
;+
;NAME:
;  GETNCDF
;
;PURPOSE: TO PERMIT EASY ACCESS OF HYPERSLABS OF GRIDDED
;			netCDF STORED IN LON,LAT, TIME or
;			LON,LAT,LEVEL (DEPTH),TIME FORMAT
;
;CATEGORY:
;	 data retrieval
;
;CALLING SEQUENCE
;
;	 GETNCDF,data,time,file,var_name,start_day,end_day,$
;		 options
;INPUTS:
;	 FILE: ID number from a call to NCDF_OPEN
;	 VAR_NAME: string name of netCDF variable to retrieve
;	 START_DAY: julian day to begin access
;	 END_DAY: julian day to end retrieval
;ex:
;	file=NCDF_OPEN('/Datasets/ncep.reanalysis.dailyavgs/surface/slp.1996.nc')
;	var_name='slp'
;	start_day=JULDAY(1,1,1996)
;  	end_day=JULDAY(1,2,1996)
;   GETNCDF,data,time,file,var_name,start_day,end_day
;		 
;		 would retrieve 2 maps for daily data
; 	
;	information on NCEP reanalysis data stored in netCDF
;	format at CDC can be found at
;	http://www.cdc.noaa.gov/cdc/reanalysis/reanalysis.shtml
;
;OUTPUTS:
;  DATA returns the hyperslab of requested netCDF data
;  TIME returns decimal true julian days 
;
;OPTIONAL KEYWORD INPUT PARAMETERS
;  NORTH: northern most latitude to be retrieved (in degrees)
;  SOUTH: southern most latitude to be retrieved (in degrees)
;  WEST: western most longitude to be retrieved  (in degrees east)    
;  EAST: eastern most longitude to be retrieved  (in degrees east)
;	 if none are specified then 0,357.5 and 90N to 90S is retrieved
;  HIGHLEVEL: largest vertical level in the dataset to be retrieved
;		 (if it has multiple levels)
;		The netCDF file must have a named variable
;			 of "level" or "depth" or "dept" or "layer"
;  LOWLEVEL: smallest vertical level in the dataset retrieved
;
;MORE OPTIONAL KEYWORD PARAMETERS
;  ANNUAL: if set, then three harmonics 
;			of annual cycle are calculated
;		   the harmonic coefficients must already exist
;		   in the file being accessed
;  AVLATITUDE: if set, averages over 
; 			the NORTH to SOUTH specified latitude range
;  AVLONGITUDE: if set, averages over the WEST to EAST specified longitude
;			   range
;  INC_LAT : interval to retrieve in latitude in units of ncdf file's latitude
;			 Ex: a file has interval 2.5 degrees. INC_LAT = 2 will
;			     retrieve every 5 degrees. Only set if need 
;			     different from 1
;  INC_LON ; interval of longitude. Only set if need different from 1
;  INC_TIME ; interval of time in units of delta_t.  Only set if need different from 1
;			inc_time will only retrieve the integer multiple of records
;			less than the total record length becuase of a bug
;			in IDL's impletation of STRIDE
;  INC_LEV; interval of levels. Only set if need different from 1
;  LEVEL: if set to levels present in the data, only the
; 	requested levels are returned. If set to a named variable 
;	returns those levels requested via highlevel,lowlevel,and inc_lev
;  cmean: if set, the climatological cmean is returned to a 
;		 named variable
;  NOcmean: if set, then the cmean is not used in calculating the climatology
;  ONLY_CLIM: if set, only the climatology is returned
;  ORIG_TIME: if set to an array, the climatology is computed
;		using these time values
;
;  TIME_UNITS: if set to a character string, the climatology is computed
;		for the orig_time array assuming these units
;  DAY_CONV: if set to original time array units conversion to 1 day
;		the climatology is computed relative to this conversion
;		
;  SPATIAL_AVERAGE: if set, spatially averages each time 
;					over latitude and longitude
;
;OPTIONAL OUTPUT PARAMETERS
;  CLIMATE: returns the 3 harmonic fit to the annual cycle
;			if ANNUAL  is set
;		   will be the same size as data. Setting ANNUAL returns
;		   returns only the annual cycle.
;	     If orig_time,time_units,day_conv are set
;		then the climate is returned relative to those time parameters 
;  LATITUDE: returns the requested latitudes
;  LONGITUDE: returns the requested longitudes
;  LEVEL: if set to levels present in the data, only the
; 	requested levels are returned. If set to a named variable 
;	returns those levels requested via highlevel,lowlevel,and inc_lev
;  MISSING_VALUE: if set to a named variable, returns the file's
;		missing_value
;  CMEAN: if set, the climatological cmean is returned to a 
;		 named variable
;  ORIG_TIME: if set to a named variable, returns the
; 		original time array as stored in the file
;
;  TIME_UNITS: if set to a named variable, returns the
;		original time array units as a string
;  DAY_CONV: if set to a named variable, returns the
;		original time array units conversion to 1 day

;  UNITS:  returns a string with the units of the data
;  TITLE:  returns a string with a short description of the data
;  VERBOSE: provides output of how and what is being retrieved
;
; NOTE: Currently fails if the dimension and variable names
;	 were written as capital letters
;
;MODIFICATION HISTORY
;	written by G. Compo on May 28, 1997
;	added inc_lat,inc_lon,inc_time 6/19/97 G. Compo
;	added cmean		8/11/97 G. Compo
;	added level keywords  2/20/98 G. Compo
;	added dept as a level option 7/15/98 G. Compo
;	changed level keyword to return non-uniform levels 11/19/98 G. Compo
;	check time units string to find offset of time array 5/13/99 G. Compo
;	removed diurnal keyword 7/19/99 G. Compo
;	annual keyword works with multi-level data 7/19/99 G. Compo
;	fixed bugs in local calculation of climatology 9/1/99 G. Compo
;	changed messages with multilevel retrieval 7/7/00 G. Compo
;	fixed bug in handling of base Julian Day 11/8/00 G. Compo
;	corrected handling of packed missing vals 3/29/01 G. Compo
;	allow for Non-COARDS time units				4/23/01 G. Compo
;	packed missing vals can be LONG even if packing is SHORT 4/26/01 G. Compo
;	allowed for a 0 year base-year, even though no year 0 in civil 
;		calendar 7/19/01 G. Compo
;	
;
; If this program is used for research purposes, please acknowledge:
; Gilbert P. Compo, NOAA/CIRES Climate Diagnostics Center,
; University of Colorado, Boulder, Colorado.
;
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
;Gil Compo                               NOAA/CIRES Climate Diagnostics Center
;Research Scientist                      325 Broadway   R/CDC1
;University of Colorado                  Boulder, CO 80303-3328
;Cooperative Institute for Research      gpc@cdc.noaa.gov
;in Environmental Sciences               (303) 497-6115,Fax:(303) 497-6449
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
;
;-
;
;///////////////////////////////////////////////////////////
PRO GETNCDF,data,time,file,var_name,start_day,end_day,$
		NORTH=north,SOUTH=south,WEST=west,EAST=east,$
		LEVEL=level,LOWLEVEL=lowlevel,HIGHLEVEL=highlevel,$
		SIG995=sig995,SFC=sfc,$
		CLIMATE=climate,ANNUAL=annual,$
		CMEAN=cmean,NOCMEAN=nocmean,$
		ORIG_TIME=orig_time,TIME_UNITS=time_units,$
		DAY_CONV=day_conv,$
		UNITS=units,TITLE=title,LATITUDE=latitude,LONGITUDE=longitude,$
		AVLATITUDE=avlatitude,AVLONGITUDE=avlongitude,$
		MISSING_VALUE=missing_value,VERBOSE=verbose,ONLY_CLIM=only_clim,$
		INC_TIME=inc_time,INC_LAT=inc_lat,INC_LON=inc_lon,$
		INC_LEV=inc_lev,TEST_OUT=test_out
		
	ON_ERROR,2 ; returns to the main level

	; check to see if levels have been specified
	IF N_ELEMENTS(level) NE 0 THEN levelspec=1
	
	; check to see if all of the time conversion
	; keywords have been set
	IF (N_elements(orig_time) ne 0) and (N_elements(day_conv) ne 0) and $
		(N_elements(time_units) ne 0) THEN BEGIN
			use_orig_time=1
	ENDIF ELSE BEGIN
		check=	(N_elements(orig_time) ne 0) +$
			 (N_elements(day_conv) ne 0) +$
			 (N_elements(time_units) ne 0)
		IF check ne 0 and check lt 3 and keyword_set(climate) THEN BEGIN
			MESSAGE,'must set all three of orig_time,day_conv, and time units',/INFO
			MESSAGE,'to use for climatology reconstruction with the orig_time passed in array'
		ENDIF
	ENDELSE
			
	;print,'THIS IS THE VAR_NAME ',var_name	
	; get the packing 
	
	gquest=NCDF_INQUIRE(file)
	never_found=1
	no_missing=1
	FOR nr=0,gquest.nvars-1 DO BEGIN
		quest=NCDF_VARINQ(file,nr)
		IF quest.name eq var_name THEN BEGIN
		  never_found=0
		  natt=quest.natts
		  FOR na=0,natt-1 DO BEGIN
			  attname=NCDF_ATTNAME(file,var_name,na)
			  CASE (attname) OF 
				  'scale_factor':BEGIN						
						  NCDF_ATTGET,file,var_name,'scale_factor',scale_factor
					  END
				   'add_offset':BEGIN						
						  NCDF_ATTGET,file,var_name,'add_offset',add_offset
					  END

				  ; if a missing value exists get its type
				  'missing_value':BEGIN
						  res=NCDF_ATTINQ(file,var_name,'missing_value')
						  miss_type=res.datatype
						  res=NCDF_VARINQ(file,var_name)
						  pack_type=res.datatype
						  NCDF_ATTGET,file,var_name,'missing_value',missing_value
						  no_missing=0
					  END
				  ELSE: BEGIN
						  IF N_ELEMENTS(scale_factor) eq 0 THEN $
							  scale_factor=1.0
						  IF N_ELEMENTS(add_offset) eq 0 THEN $
							  add_offset=0.0
					  END
			  ENDCASE
		  ENDFOR
		ENDIF
	ENDFOR
	IF never_found THEN BEGIN
		MESSAGE,'NO variable by the name '+var_name+' in this file',/INFO
		
		MESSAGE,'Here are the variable names in this file ',/INFO
		FOR nu=0,gquest.Nvars-1 DO BEGIN
			aquest=NCDF_VARINQ(file,nu)
			MESSAGE,aquest.name,/INFO
		ENDFOR
		MESSAGE,'please try one of the above names'
	ENDIF
		
	
	; if haven't found a missing_value, look under global
	IF no_missing THEN BEGIN
	
		;look for a global missing value
		FOR nn=0,gquest.ngatts-1 DO BEGIN
			attname=NCDF_ATTNAME(file,/GLOBAL,nn)
			CASE (attname) OF 
				; if a missing value exists get its type
				'missing_value':BEGIN
						res=NCDF_ATTINQ(file,/GLOBAL,'missing_value')
						miss_type=res.datatype
						res=NCDF_VARINQ(file,var_name)
						pack_type=res.datatype
						NCDF_ATTGET,file,/GLOBAL,'missing_value',missing_value
						END
				ELSE: BEGIN
					END
			ENDCASE
		ENDFOR	
	
	ENDIF
	;print,add_offset,scale_factor,missing_value
	
	; check to see if latitude, longitude, level, or time increments are set
	IF NOT(KEYWORD_SET(inc_lat)) THEN inc_lat=1
	IF NOT(KEYWORD_SET(inc_lon)) THEN inc_lon=1
	IF NOT(KEYWORD_SET(inc_time)) THEN inc_time=1
	IF NOT(KEYWORD_SET(inc_lev)) THEN inc_lev=1
	; get the ancillary data
	res=NCDF_INQUIRE(file)
	nvar=res.nvars
	FOR nn=0,nvar-1 DO BEGIN
		quest=NCDF_VARINQ(file,nn)
		IF quest.name eq var_name THEN BEGIN
			natt=quest.natts
			ndims=quest.ndims
			dim_arr=quest.dim
			FOR na=0,natt-1 DO BEGIN
				attname=NCDF_ATTNAME(file,var_name,na)
				CASE (attname) OF 
					'long_name':BEGIN
								NCDF_ATTGET,file,var_name,'long_name',title 
								title=STRING(title)
								END
					'title':BEGIN
								NCDF_ATTGET,file,var_name,'title',title 
								title=STRING(title)
								END
					'units':BEGIN
								NCDF_ATTGET,file,var_name,'units',units
								units=STRING(units)
								END
					ELSE: BEGIN
								IF N_ELEMENTS(title) eq 0 THEN title=''
								IF N_ELEMENTS(units) eq 0 THEN units=''
								END

				ENDCASE
			ENDFOR

		ENDIF 
	ENDFOR
	
	; construct the time arrays. time is in days, htime is in delta_t since starting from (1,1,1)
	; htime is needed to construct the anomalies if annual_coef exists
	NCDF_DIMINQ,file,'time',name,timsiz
	NCDF_VARGET,file,'time',time0,count=[1]
	; units are in hours for some and days for others, sometimes MONTHS!
	NCDF_ATTGET,file,'time','units',tunits
	tst=STRMID(tunits,0,1) ; find out what the time units are
	IF (N_ELEMENTS(time_units) EQ 0) THEN BEGIN
		time_units=STRING(tunits)
	ENDIF ELSE BEGIN
		IF (time_units NE STRING(tunits)) AND Keyword_set(climate) THEN BEGIN
			MESSAGE,'time_units must be the same as the time units of the coefficients',/INFO
			PRINT,'time_units passed in= ',time_units
			PRINT,'time units in this file= ',STRING(tunits)
			MESSAGE,'please make these the same or unset keyword time_units'
		ENDIF
	ENDELSE
	CASE tst OF
		'h':dt=24D0
		'd':dt=1D0
		's':dt=86400D0
		'T':dt=1D0 ; T stands for True Julian Day used by NOAA PMEL for TAO buoy data
		'm':dt=30D0 ; CRU uses months since some day
	ENDCASE
	IF (N_ELEMENTS(day_conv) EQ 0) THEN day_conv=dt

	; convert to days
	; try to find if there is a base_jd attribute
	; also trying to find out the time spacing between obs 
	; not everyone follows COARDS conventions so really looking around here
	not_found_delta_t=1
	not_found_base_jd=1
	FOR nn=0,nvar-1 DO BEGIN
		quest=NCDF_VARINQ(file,nn)
		IF quest.name eq 'time' THEN BEGIN
			natt=quest.natts
			FOR na=0,natt-1 DO BEGIN
				attname=NCDF_ATTNAME(file,'time',na)
				CASE (attname) OF 
					'base_jd':BEGIN
							NCDF_ATTGET,file,'time','base_jd',base_jd 
							not_found_base_jd=0
						END
					'delta_t':BEGIN
						not_found_delta_t=0
						NCDF_ATTGET,file,'time','delta_t',delta_t
						READS,STRING(delta_t),dyr,dmnth,ddy,dh,dm,ds,$
							format='(I4,5(1X,I2,1X,I2))'
						END
					ELSE: BEGIN
							
							IF STRING(tunits) eq 'True Julian Day' THEN BEGIN
								base_jd=0
								not_found_base_jd=0
							ENDIF
						END
						
				ENDCASE
			ENDFOR

		ENDIF 
	ENDFOR
	;haven't found base_jd so convert from the time units
	IF not_found_base_jd THEN BEGIN
		; find the string position of 'since'
		tunits=STRING(tunits)
		since_pos=STRPOS(tunits,'since')
		since_off=5
		IF since_pos eq -1 THEN BEGIN
			; try from
			since_pos=STRPOS(tunits,'from')
			since_off=4
		ENDIF

		IF since_pos ne -1 THEN BEGIN
			; try to find two dashes and a colon 
			; or two slashes and a colon to
			; construct the base time
			tunits_len=STRLEN(tunits)
			sub_string=STRMID(tunits,since_pos+since_off,tunits_len-1)
			; remove leading blanks
			sub_string=STRTRIM(sub_string,1)
			sub_len=STRLEN(sub_string)
			;PRINT,sub_string
			; find first dash
			dash1=STRPOS(sub_string,'-')
			IF dash1 EQ -1 THEN BEGIN
				; try for slashes
				dash1=STRPOS(sub_string,'/') 
				dash2=RSTRPOS(sub_string,'/')
			ENDIF ELSE BEGIN
				dash2=RSTRPOS(sub_string,'-')
			ENDELSE
			IF (dash1 NE -1) AND (dash2 NE -1) THEN BEGIN
				IF dash1 LT 4 THEN BEGIN
					month=FLOAT(STRMID(sub_string,0,dash1-0))
					;PRINT,month
					day=FLOAT(STRMID(sub_string,dash1+1,dash2-dash1-1))
					;find the final space
					end_ymd_mark=STRPOS(sub_string,' ')
					IF end_ymd_mark ne -1 THEN $
						year=FLOAT(STRMID(sub_string,dash2+1,end_ymd_mark-dash2-1)) ELSE $
						year=FLOAT(STRMID(sub_string,dash2+1,sub_len-dash2-1))
					;PRINT,month,day,year
				ENDIF ELSE BEGIN
					year=FLOAT(STRMID(sub_string,0,dash1-0))
					month=FLOAT(STRMID(sub_string,dash1+1,dash2-dash1))
					end_ymd_mark=STRPOS(sub_string,' ')
					IF end_ymd_mark ne -1 THEN $
						day=FLOAT(STRMID(sub_string,dash2+1,end_ymd_mark-dash2-1)) ELSE $
						day=FLOAT(STRMID(sub_string,dash2+1,sub_len-dash2-1))
					; assume that the year was placed first
				ENDELSE	

				; If end_ymd_mark of a space does not exist then no need to look for hour,minute,second offset
				hour=0.0d0
				minute=0.0d0
				second=0.0d0
				IF end_ymd_mark LT sub_len AND end_ymd_mark ne -1 THEN BEGIN
					colon1=STRPOS(sub_string,':')
					colon2=RSTRPOS(sub_string,':')
					hour=FLOAT(STRMID(sub_string,colon1-2,2)) ; hours can only be two digits
					minute=FLOAT(STRMID(sub_string,colon1+1,2)) ; minutes can only be two digits
					second=FLOAT(STRMID(sub_string,colon2+1,sub_len-colon2))
					;PRINT,hour,minute,second
				ENDIF
				; there is no year 0 in the civil calendar
				; but NCAR has indicated a year 0 in several model runs,
				; so be read for it
				IF year LT 0.9 THEN year=1
				base_jd=JULDAY(month,day,year)+hour/24.0d0+minute/(24.0d0*60.0d0)+second/(86400.0d0) 

				not_found_base_jd=0
			ENDIF 

		ENDIF 
	ENDIF
	
	
	; haven't found a delta_t, so need to construct what the time spacing is between obs
	IF not_found_delta_t THEN BEGIN
		; set the deltas all to zero
		dyr=0.0 & dmnth=0.0 &  ddy=0.0 & dh=0.0 & dm=0.0 & ds=0.0
		
		; look for a time2
		FOR nn=0,nvar-1 DO BEGIN
		  quest=NCDF_VARINQ(file,nn)
		  IF quest.name eq 'time2' THEN BEGIN
			  natt=quest.natts
			  FOR na=0,natt-1 DO BEGIN
				  attname=NCDF_ATTNAME(file,'time2',na)
				  CASE (attname) OF 
					  'units':BEGIN ; right now just checking for TAO conventions
						  NCDF_ATTGET,file,quest.name,attname,t2units
						  t2units=STRING(t2units)
						  strloc=STRPOS(t2units,'msec')
					    	  IF strloc ne -1 THEN BEGIN
						     NCDF_VARGET,file,'time2',time2,count=2
						     time2=time2/1000. ; now in seconds
						     ds=time2(1)-time2(0)
						     IF ds eq 0 THEN BEGIN
							   ; assume daily
							   NCDF_VARGET,file,'time',utime,count=2
							   ddy=utime(1)-utime(0)
						     ENDIF
						  ENDIF
						  END
					    ELSE:BEGIN
					    	 END

				  ENDCASE
			  ENDFOR

		  ENDIF 
		ENDFOR	
		
		
	ENDIF 						
	delta_t=(dyr*365.25*DT+dmnth*30.0*DT+ddy*DT+dh*DT/24.+dm*DT/(60.*24.0)+ds*DT/(24.*60.*60.))
	htime=time0+dindgen(timsiz)*delta_t
	
	; this should make the time units in days
	 IF var_name eq 'sst' THEN BEGIN 
	 	time0=time0+base_jd  ; need this for Reynold's SST data
	 ENDIF ELSE BEGIN
	 	IF not_found_base_jd THEN BEGIN
			time0=time0/dt 
			MESSAGE,'COULD NOT FIGURE OUT THE TIME UNITS AND SPACING',/INFO
			MESSAGE,'SO I WILL RETURN ALL OF THE TIME DATA',/INFO
			MESSAGE,'UNLESS start_day and end_day correspond to days in the',/INFO
			MESSAGE,'increments of the time variable',/INFO
		ENDIF ELSE BEGIN
			time0=time0/dt+base_jd
		ENDELSE
	 ENDELSE
	
	
	; want the time array to be in days no matter what the original units
	delta_t=(dyr*365.25+dmnth*30.0+ddy+dh/24.+dm/(24.*60.)+ds/(24.*60.*60.))
	time=dindgen(timsiz)*delta_t
	time=time+time0

	
	
	; if have weekly or monthly data, then need the time arrays from the
	; original file instead because julian days will not be evenly spaced (potentially)
	IF delta_t ge 7.0 OR not_found_delta_t THEN BEGIN
		IF keyword_set(verbose) THEN BEGIN
			MESSAGE,'my delta_t is greater than 7.0 ',/INFO
			MESSAGE,'OR I could not find a delta_t',/INFO
			MESSAGE,'SO I am doing my best to get the correct time',/INFO
			MESSAGE,'I recommend checking the results carefully',/INFO
		ENDIF
		NCDF_VARGET,file,'time',time
		htime=time
		time=DOUBLE(time)
		time=time/dt+base_jd ; puts in julian days since base_jd
		
		; perhaps the data is monthly and so days are not equally spaced!
		; and no delta_t was given
		IF STRMID(tunits,0,3) eq 'mon' and not_found_base_jd eq 0 THEN BEGIN
			IF keyword_set(verbose) THEN MESSAGE,'I think I have monthly data',/INFO
			; I assume I was able to figure out the Julian day of the first time
			nmonths=timsiz
			CALDAT,base_jd,umonth,udy,uyr

			FOR jt=0,nmonths-1 DO BEGIN
				time[jt]=JULDAY(umonth,1,uyr)
				umonth=umonth+1
				IF umonth GT 12 THEN BEGIN
					umonth=1
					uyr=uyr+1
				ENDIF
			ENDFOR;jt
			
		ENDIF
	ENDIF

	IF not_found_base_jd THEN BEGIN
		start_day=MAX([MIN(time),start_day])
		end_day=MIN([MAX(time),end_day])
	ENDIF
	
	; starting location of hyperslab in time
	ttot=where(time ge start_day and time le end_day)
	tint=LONG(N_ELEMENTS(ttot))
	IF inc_time ge tint THEN inc_time=1
	
	ttot=ttot(LINDGEN(tint/inc_time)*inc_time)
	t=MIN(ttot)
	tint=LONG(N_ELEMENTS(ttot)) 
	htime=htime(ttot)
	time=time(ttot)
	
	; construct the levels requested
	
	;find out the names of the dimensions in the variable that we are pursuing
	
	; first find out if levels or depth variables
	yeslev=0
	FOR nn=0,nvar-1 DO BEGIN
		quest=NCDF_VARINQ(file,nn)
		IF quest.name eq 'level' or quest.name eq 'lev' or $
		quest.name eq 'layer' or  $
		(STRPOS(quest.name,'dep') NE -1) THEN BEGIN
		  ; check against the dimensions of the var_name to see if our var
		  ; has this level variable as a coordinate dimension
		  	FOR nd=0,ndims-1 DO BEGIN
				dimid=dim_arr(nd)
				NCDF_DIMINQ,file,dimid,dimname,dimsize
				IF quest.name eq dimname THEN BEGIN
					NCDF_VARGET,file,quest.name,alevels
					yeslev=1
				ENDIF	
			ENDFOR		
		ENDIF
	ENDFOR
	IF yeslev THEN BEGIN
		IF (N_ELEMENTS(highlevel) eq 0) THEN highlevel=MAX(alevels)
		IF (N_ELEMENTS(lowlevel) eq 0) THEN lowlevel=MIN(alevels)
		IF (lowlevel gt highlevel) THEN $
			MESSAGE,'LOWLEVEL must be less than or equal to HIGHLEVEL'
			
		IF (N_ELEMENTS(levelspec) eq 0) THEN BEGIN
			uselevel=WHERE(alevels le highlevel and alevels ge lowlevel) 
		ENDIF ELSE BEGIN
			highlevel=MAX(level) & lowlevel=MIN(level)
			uselevel=WHERE(alevels le highlevel and alevels ge lowlevel) 
			wantlevel=REMIX(level,alevels(uselevel))
			IF KEYWORD_SET(verbose) AND N_Elements(wantlevel) GT 1 THEN BEGIN
				MESSAGE,'all levels in the range will be retrieved and '+$
				'the requested subset will be returned',/INFO
				MESSAGE,'To retrieve fewer levels, ',/INFO
				MESSAGE,'employ keywords highlevel,lowlevel, and inc_level ',/INFO
			ENDIF
			
		ENDELSE
		IF uselevel(0) eq -1 THEN BEGIN
			IF (N_ELEMENTS(levelspec) eq 0) THEN BEGIN
				MESSAGE,'LOWLEVEL='+STRTRIM(STRING(lowlevel),1),/info
				MESSAGE,'HIGHLEVEL='+STRTRIM(STRING(highlevel),1),/info
				MESSAGE,'ACTUAL RANGE OF LEVELS '+$
					STRTRIM(STRING(MAX(alevels)),1)+','+$
					STRTRIM(STRING(MIN(alevels)),1),/info
				MESSAGE,'LOWLEVEL and HIGHLEVEL must be in range of available levels'
			ENDIF ELSE BEGIN
				MESSAGE,'ACTUAL LEVELS ',/info
				PRINT,	STRTRIM(STRING(alevels),1)
				MESSAGE,'YOUR REQUESTED LEVELS ',/info
				PRINT,STRTRIM(STRING(level),1)
				MESSAGE,'REQUESTED LEVELS MUST BE A SUBSET OF ACTUAL LEVELS'
			ENDELSE
		ENDIF
		nlev=N_ELEMENTS(uselevel)
		IF inc_lev ge nlev THEN inc_lev=1 ; cannot stride more than minimum levels requested
		IF (N_ELEMENTS(levelspec) eq 0) THEN inc_lev=1
		uselevel=uselevel(LINDGEN(nlev/inc_lev)*inc_lev)
		 
		nlev=N_ELEMENTS(uselevel)
		stlev=MIN(uselevel)
		hold_level=alevels(uselevel)
	ENDIF ELSE BEGIN
		IF ((N_ELEMENTS(highlevel) NE 0) OR (N_ELEMENTS(lowlevel) NE 0) OR $
			(N_ELEMENTS(levelspec)  NE 0)) THEN $
			MESSAGE,'LEVEL REQUESTS IGNORED BECAUSE NO LEVEL-type dimension',/info
	ENDELSE
	
	; construct the latitudes and longitudes desired
	NCDF_VARGET,file,'lat',lat
	NCDF_VARGET,file,'lon',lon
	IF (N_ELEMENTS(north) eq 0) THEN north=max(lat)
	IF (N_ELEMENTS(south) eq 0) THEN south=min(lat)
	IF (N_ELEMENTS(west) eq 0) THEN west=min(lon)
	IF (N_ELEMENTS(east) eq 0) THEN east=max(lon)
	IF (west gt east) THEN MESSAGE,'WEST must be less or equal to EAST'
	uselat=where(lat le north and lat ge south) 
	nlat=N_ELEMENTS(uselat)
	IF inc_lat ge nlat THEN inc_lat=1 ;  cannot stride more than the minimum distance 
									  ; requested
	uselat=uselat(LINDGEN(nlat/inc_lat)*inc_lat)
	latitude=lat(uselat)
	
	uselon=where(lon ge west and lon le east) 
	nlon=N_ELEMENTS(uselon)
	IF inc_lon ge nlon THEN inc_lon=1 ;  cannot stride more than the minimum distance 
									  ; requested
	uselon=uselon(LINDGEN(nlon/inc_lon)*inc_lon)
	longitude=lon(uselon)
	nlat=N_ELEMENTS(uselat)
	nlon=N_ELEMENTS(uselon)
	stlat=MIN(uselat)
	stlon=MIN(uselon)
	
	;SET the offsets
	IF yeslev THEN BEGIN
		count=[nlon,nlat,nlev,tint]
		offset=[stlon,stlat,stlev,t]
		stride=[inc_lon,inc_lat,inc_lev,inc_time]
		;PRINT,'these are the counts, offset, and stride being used'
		;PRINT,count
		;PRINT,offset
		;PRINT,stride
	ENDIF ELSE BEGIN
		count=[nlon,nlat,tint]
		offset=[stlon,stlat,t]
		stride=[inc_lon,inc_lat,inc_time]
	ENDELSE
	
	; CHECK TO SEE THAT THE DATA AMOUNT REQUESTED IS NOT TOO LARGE
	amount_req=(FLOAT(nlon)/inc_lon)*(FLOAT(nlat)/inc_lat)*(FLOAT(tint)/inc_time)
	amount_req=amount_req*8L
	IF amount_req gt 1.0E+9 and NOT(KEYWORD_SET(only_clim)) THEN BEGIN
			MESSAGE,'YO! YOU HAVE REQUESTED '+$
					ROUND_ANY(amount_req)+' MB OF MEMORY',/INFO
			MESSAGE,'THIS IS EXCESSIVE. PLEASE REDUCE YOUR REQUEST'
	ENDIF
		
	; get the data
	IF keyword_set(test_out) THEN STOP

	IF KEYWORD_SET(verbose) and NOT(KEYWORD_SET(only_clim)) THEN BEGIN
		print,' '
		print,'I am getting the data of dimensions ',count
		print,' beginning at position longitude =',lon(offset(0)), ' latitude=',lat(offset(1))
		print,'STRIDE =',stride
		PRINT,'COUNT =',count
		PRINT,'OFFSET =',offset
		print,'here is the info on the var_name, count,offset, and stride being used'
		hlp,file,var_name,count,offset,stride
		IF tint ge 1000 THEN $
			print,'If your request is a long time array, this could really take awhile. STAY CALM :-]'
	ENDIF
	IF NOT(KEYWORD_SET(only_clim)) THEN BEGIN
		NCDF_VARGET,file,var_name,data,count=count,offset=offset,stride=stride
	ENDIF ELSE BEGIN
		data=0.
	ENDELSE
	IF KEYWORD_SET(verbose) THEN print,'I have gotten the data and it will be to you soon :)'
	IF (KEYWORD_SET(climate)) THEN BEGIN
	
	; Check to see if annual cycle coefficients exist
		res=NCDF_INQUIRE(file)
		nvar=res.nvars
		noharmfound=1
		FOR nn=0,nvar-1 DO BEGIN
			quest=NCDF_VARINQ(file,nn)
			
			IF quest.name eq 'annual_coef' THEN BEGIN ; only annual_coef variable for 
								;multiple level data
				;print,'GOT the annual coefs'
				noharmfound=0
				NCDF_DIMINQ,file,'harm_coef',name,ncoef
				IF yeslev THEN BEGIN
				NCDF_VARGET,file,'annual_coef',annual_coef,$
					stride=[inc_lon,inc_lat,inc_lev,1],$
					count=[nlon,nlat,nlev,ncoef],$
					offset=[stlon,stlat,stlev,0]
				ENDIF ELSE BEGIN
				NCDF_VARGET,file,'annual_coef',annual_coef,$
					stride=[inc_lon,inc_lat,1],count=[nlon,nlat,ncoef],offset=[stlon,stlat,0]
				ENDELSE
			ENDIF ELSE BEGIN
				IF N_ELEMENTS(ncoef) eq 0 THEN ncoef=7
				IF (N_elements(annual_coef) EQ 0) THEN $
					IF yeslev THEN annual_coef=fltarr(nlon,nlat,nlev,ncoef) ELSE $
					annual_coef=fltarr(nlon,nlat,ncoef)
			ENDELSE
			IF quest.name eq 'harm_coef_'+var_name THEN BEGIN
				noharmfound=0
				NCDF_DIMINQ,file,'harm_coef',name,ncoef
				IF yeslev THEN BEGIN
				NCDF_VARGET,file,'harm_coef_'+var_name,annual_coef,$
					stride=[inc_lon,inc_lat,inc_lev,1],$
					count=[nlon,nlat,nlev,ncoef],$
					offset=[stlon,stlat,stlev,0]
				ENDIF ELSE BEGIN
				NCDF_VARGET,file,'harm_coef_'+var_name,annual_coef,$
					stride=[inc_lon,inc_lat,1],count=[nlon,nlat,ncoef],offset=[stlon,stlat,0]
				ENDELSE
			ENDIF ELSE BEGIN
				IF (n_elements(annual_coef) EQ 0) THEN BEGIN
					print,'screwing with the annual coefs'
					IF N_ELEMENTS(ncoef) eq 0 THEN ncoef=7
					IF yeslev THEN annual_coef=fltarr(nlon,nlat,nlev,ncoef) ELSE $ 
					annual_coef=fltarr(nlon,nlat,ncoef)
				ENDIF
			ENDELSE
		ENDFOR
		
		IF (KEYWORD_SET(cmean)) THEN BEGIN
				IF yeslev THEN cmean=REFORM(annual_coef(*,*,*,0)) ELSE $
					cmean=REFORM(annual_coef(*,*,0))
		ENDIF
		IF (KEYWORD_SET(nocmean)) THEN BEGIN
			IF yeslev THEN BEGIN
				annual_coef(*,*,*,0)=0.0
			ENDIF ELSE BEGIN
				annual_coef(*,*,0)=0.0
			ENDELSE
		ENDIF
		
		;annual_coef=annual_coef(uselon,*,*) & annual_coef=annual_coef(*,uselat,*)

		
		;hlp,annual_coef,diurnal_coef
		IF KEYWORD_SET(verbose) THEN BEGIN
			print,'computing the harmonic cycles'
			print,'using '+ROUND_ANY(ncoef)+' coefficients '
		ENDIF
		IF noharmfound THEN BEGIN
			IF KEYWORD_SET(verbose) THEN BEGIN
				MESSAGE,'no harmonic coefficients were ever found',/INFO
			ENDIF
		ENDIF ELSE BEGIN
			IF (N_ELEMENTS(orig_time) eq 0) THEN BEGIN	
				CLIMSER,htime,dt,annual_coef,(ncoef-1),climate
			ENDIF ELSE BEGIN
				IF Keyword_set(verbose) THEN $
					MESSAGE,'using the time array passed in with orig_time',/INFO
				CLIMSER,orig_time,day_conv,annual_coef,(ncoef-1),climate
			ENDELSE
		ENDELSE
	ENDIF;keyword set climate
	; check to see if any missing values present in the dataset

	IF N_ELEMENTS(missing_value) NE 0 THEN BEGIN
	 IF miss_type eq pack_type THEN BEGIN
		missloc=where(ABS(float(data)-float(missing_value)) LT 1.0E-6)
		data=TEMPORARY(data)*scale_factor+add_offset
		
	 ENDIF ELSE BEGIN
	 	; now check and see if pack type is short and
		; the missing is a LONG
		IF (miss_type eq 'LONG' AND pack_type EQ 'SHORT') OR $
			(miss_type eq 'SHORT' AND pack_type EQ 'LONG') THEN BEGIN
			
			missloc=where(ABS(float(data)-float(missing_value)) LT 1.0E-6)
			data=TEMPORARY(data)*scale_factor+add_offset
		ENDIF ELSE BEGIN
		
	; means that missing_value must be the same type as the unpacked data
		data=TEMPORARY(data)*scale_factor+add_offset
		missloc=where(ABS(data-missing_value) LT 1E-6)
		ENDELSE
		
	 ENDELSE
	 IF (missloc(0) ne -1) THEN data(missloc)=missing_value
	 
	ENDIF ELSE BEGIN
		data=TEMPORARY(data)*scale_factor+add_offset
	ENDELSE

	;hlp,data
	
	IF ((N_ELEMENTS(level) NE 0) and yeslev) THEN BEGIN
	
		data=data(*,*,wantlevel,*)
		level=alevels(uselevel(wantlevel))

		IF (KEYWORD_SET(climate)) THEN BEGIN

		   IF noharmfound THEN BEGIN ; if want a climatology, calculate it
					   ; for the case of no coefficients
			   IF KEYWORD_SET(verbose) THEN BEGIN
				   MESSAGE,'no harmonic coefficients were ever found',/INFO
				   MESSAGE,'But climatology keyword is set',/INFO
				   MESSAGE,'so computing the 3 harmonic annual cycle',/INFO
			   ENDIF

			   IF yeslev THEN BEGIN
				   cmean=TOTAL(data,4)/tint
				   climate=data
				   FOR i=0,nlon DO BEGIN
				   FOR j=0,nlat DO BEGIN
				   FOR k=0,N_ELEMENTS(wantlevel)-1 DO BEGIN
					   climate(i,j,k,*)=HARMFIT(time,data(i,j,k,*)-cmean(i,j,k),$
								   365.25D0,3)
				   ENDFOR;k
				   ENDFOR;j
				   ENDFOR;i

			   ENDIF ELSE BEGIN
				   cmean=TOTAL(data,3)/tint
				   climate=data
				   FOR i=0,nlon DO BEGIN
				   FOR j=0,nlat DO BEGIN
					   climate(i,j,*)=HARMFIT(time,data(i,j,*)-cmean(i,j),$
								   365.25D0,3)
				   ENDFOR;
				   ENDFOR;

			   ENDELSE
		   ENDIF ELSE BEGIN
			   climate=climate(*,*,wantlevel,*)
		   ENDELSE
		ENDIF
	ENDIF ELSE BEGIN
		IF yeslev THEN level=alevels(uselevel)
		
		IF (KEYWORD_SET(climate)) THEN BEGIN

		   IF noharmfound THEN BEGIN ; if want a climatology, calculate it
					   ; for the case of no coefficients
			   IF KEYWORD_SET(verbose) THEN BEGIN
				   MESSAGE,'no harmonic coefficients were ever found',/INFO
				   MESSAGE,'But climatology keyword is set',/INFO
				   MESSAGE,'so computing the 3 harmonic annual cycle',/INFO
			   ENDIF

			   IF yeslev THEN BEGIN
				   cmean=TOTAL(data,4)/tint
				   climate=data
				   FOR i=0,nlon-1 DO BEGIN
				   FOR j=0,nlat-1 DO BEGIN
				   FOR k=0,N_ELEMENTS(wantlevel)-1 DO BEGIN
					   climate(i,j,k,*)=HARMFIT(time,data(i,j,k,*)-cmean(i,j,k),$
								   365.25D0,3)
								   
				   ENDFOR;k
				   ENDFOR;j
				   ENDFOR;i
				   IF NOT(keyword_set(nocmean)) THEN climate=TEMPORARY(climate)+$
				   	REBIN(REFORM(cmean,nlon,nlat,nlevel,1),nlon,nlat,nlevel,tint)

			   ENDIF ELSE BEGIN
				   cmean=TOTAL(data,3)/tint
				   climate=data
				   FOR i=0,nlon-1 DO BEGIN
				   FOR j=0,nlat-1 DO BEGIN
					   climate(i,j,*)=HARMFIT(time,data(i,j,*)-cmean(i,j),$
								   365.25D0,3)
								  
				   ENDFOR;
				   ENDFOR;
				   IF NOT(keyword_set(nocmean)) THEN climate=TEMPORARY(climate)+$
				   	REBIN(REFORM(cmean,nlon,nlat,1),nlon,nlat,tint)

			   ENDELSE
		   ENDIF 
		ENDIF
		
	ENDELSE
	 
	IF (KEYWORD_SET(avlatitude) or KEYWORD_SET(avlongitude)) THEN BEGIN
		IF (KEYWORD_SET(avlatitude) and KEYWORD_SET(avlongitude)) THEN BEGIN
			data=AREA_AV(data)
			IF (KEYWORD_SET(climate)) THEN climate=AREA_AV(climate)
		ENDIF  ELSE BEGIN
			IF (KEYWORD_SET(avlatitude)) THEN BEGIN
				print,'Averaging in latitude'
				data=TOTAL(TEMPORARY(data),2)/nlat
				IF (KEYWORD_SET(climate)) THEN climate=TOTAL(TEMPORARY(climate),2)/nlat
			ENDIF ELSE BEGIN
				IF (KEYWORD_SET(avlongitude)) THEN BEGIN
					print,'Averaging in longitude'
					data=TOTAL(TEMPORARY(data),1)/nlon
					IF (KEYWORD_SET(climate)) THEN climate=TOTAL(TEMPORARY(climate),1)/nlon
				ENDIF
			ENDELSE
		ENDELSE
	ENDIF
	
	; set the original time for returning
	IF N_ELEMENTS(orig_time) EQ 0 THEN orig_time=htime
	IF yeslev THEN level=alevels(uselevel)
	If N_ELEMENTS(levelspec) NE 0 AND yeslev THEN level=hold_level(wantlevel)

	
	RETURN
	END