extern double exp(), log();
extern float getsltlon();

/*
 * Return indices n1,n2 to xx, bracketing value x
 */
int bracket(x,xx,nx,n1,n2)
float x, *xx; 
int nx, *n1, *n2;
{
	int i;
	if (x < *xx) { *n1 = 0; *n2 = 1; return(0); }
	if (x > *(xx+nx-1)) { *n1 = nx-2; *n2 = nx-1; return(0); }
	for (i=0; i < nx; i++) {
	  if (x >= *(xx+i) && x <= *(xx+i+1)) {
	    *n1 = i; *n2 = i+1; return(0);
	  }
	}
	return(-1);
}

/*
 * Interpolate gcmin(kmx) (a field pressure column) to hts(nhts) heights, 
 * returning interpolated values in gcmout(nhts). If loght > 0, do log 
 * (exponential) interpolation, linear otherwise. Use gcmht(kmx) heights 
 * (heights in pressure column) to bracket desired heights. If desired
 * height in hts is above max gcmht or below min gcmht, return spval.
 */
int htint(gcmin,gcmht,loght,hts,nhts,kmx,gcmout,spval)
float *gcmin, *gcmht, *hts, *gcmout, spval;
int loght, nhts, kmx;
{
	int k, ih, k0,k1;
	float htmin, htmax, f1;

	htmin = 1.e36; htmax=-1.e36;
	for (k=0; k < kmx; k++) {
	  if (*(gcmht+k) < htmin) htmin = *(gcmht+k);
	  if (*(gcmht+k) > htmax) htmax = *(gcmht+k);
	}
	for (ih=0; ih < nhts; ih++) {
	  if (*(hts+ih) > htmax || *(hts+ih) < htmin) {
	    *(gcmout+ih) = spval;
	  } else {
	    bracket(*(hts+ih),gcmht,kmx,&k0,&k1);
	    if (!loght) {
              f1 = (*(hts+ih) - *(gcmht+k0)) / (*(gcmht+k1) - *(gcmht+k0));
              f1 = f1 * (*(gcmin+k1)) + (1.-f1) * (*(gcmin+k0));
	    } else {
              f1 = *(gcmin+k1) / *(gcmin+k0);
              f1 = (log((double)f1) / 
	           (*(gcmht+k1) - *(gcmht+k0))) * (*(hts+ih) - *(gcmht+k0));
              f1 = *(gcmin+k0) * exp((double)f1);
	    } /* loght */
	    *(gcmout+ih) = f1;
	  } /* if within height range */  
	} /* ih= 0, nhts-1 */
}

/*
 * Bilinearly interpolate gcmin to slat,slon, returning interpolated 
 * pressure col in gcmout. If (interpslt) then interpolate to sat 
 * slt (sslt) rather than to slon (i.e., find longitude corresponding 
 * to gcm ut and sslt, and interpolate to that longitude)
 */
int locint(slat,slon,sslt,gcmin,gcmut,interpslt,imx,kmx,jmx,gcmlat,
	   gcmlon,gcmout)
float slat,slon,sslt, *gcmin, gcmut, *gcmlat, *gcmlon, *gcmout;
int imx,kmx,jmx,interpslt;
{
	int k, lat0, lat1, lon0, lon1;
	float flon, flat, f1,f2,f3,f4, sltlon;

	bracket(slat,gcmlat,jmx,&lat0,&lat1);
	if (interpslt) {
	  sltlon = getsltlon(sslt,gcmut);
	  bracket(sltlon,gcmlon,imx,&lon0,&lon1);
	} else
	  bracket(slon,gcmlon,imx,&lon0,&lon1);
/*
 * Think of gcmin as [jmx][kmx][imx], then:
 */
	for (k=0; k < kmx; k++) {
	  if (interpslt)
	    flon = (sltlon-(*(gcmlon+lon0))) / 
	      (*(gcmlon+lon1)-(*(gcmlon+lon0)));
	  else
	    flon = (slon-(*(gcmlon+lon0))) / 
	      (*(gcmlon+lon1)-(*(gcmlon+lon0)));
	  flat = (slat-(*(gcmlat+lat0))) / (*(gcmlat+lat1)-(*(gcmlat+lat0)));
	  f1 = *(gcmin + lat0*kmx*imx + k*imx + lon0);
	  f2 = *(gcmin + lat0*kmx*imx + k*imx + lon1);
	  f3 = *(gcmin + lat1*kmx*imx + k*imx + lon0);
	  f4 = *(gcmin + lat1*kmx*imx + k*imx + lon1);
	  *(gcmout+k) = flon*(flat*f4+(1.-flat)*f2) +
	                (1.-flon)*(flat*f3+(1.-flat)*f1);
	}
	return(0);
}