/**---------------------------------------------------------------------------
 ** 
 ** vgeo.c -- 
 ** 
 ** Author          : Pierre Jaccard
 ** Created On      : 1999/08/04 14:56:17
 ** Last Modified By: Pierre Jaccard
 ** Last Modified On: 1999/08/04 17:42:30
 ** Update Count    : 6
 ** Directory       : /home/pego/pcd1/codas3c/gfi/src/apps/util/
 ** Version         : 0.0
 ** Status          : Unknown
 ** ---------------------------------------------------------------------- ** 
 ** DESCRIPTION: 
 ** 
 **    Undocumented.
 ** 
 ** ---------------------------------------------------------------------- ** 
 ** REVISIONS: 
 ** ---------------------------------------------------------------------- ** 
 ** CHANGES: 
 **------------------------------------------------------------------------**/

#include "vgeo.h"

/* 1999/08/04 Pierre Jaccard
   copy_ymdhms_time was originally located in use_db/time_io.c, but has
   disapeared with the new CODAS files from Eric.
*/
#if PROTOTYPE_ALLOWED
void copy_ymdhms_time(YMDHMS_TIME_TYPE *dest, YMDHMS_TIME_TYPE *source)
#else
void copy_ymdhms_time(dest, source)
YMDHMS_TIME_TYPE *dest, *source;
#endif
{
   dest->year   = source->year;
   dest->month  = source->month;
   dest->day    = source->day;
   dest->hour   = source->hour;
   dest->minute = source->minute;
   dest->second = source->second;
}                       /* copy_ymdhms_time() */

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Reads further options according to an optionnal keyword. The argument
	    is assumed to be a table of correspondance between a list of options and
			the value of <code>. If the table has an entry for the given <code>, the
	    corresponding options are executed.
	 
	 RETURNS: The value of <code>.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
char *vgeo_read_options(/* File descriptor to the control file */
												FILE *fp, 
												/* Pointer to the table of correspondance between
													 <code> and further options to be exectuted */
												char *arg, 
												/* Index for finding the list of options to be
													 executed */
												int code
												)
#else
char *vgeo_read_options(fp, arg, code)
FILE *fp;
char *arg;
int  code;
#endif
{

	VGEO_OPTION_TABLE_TYPE *list, *p;
	char msg[200];

	/* Check if there are specific options to execute */
	list = (VGEO_OPTION_TABLE_TYPE *) arg;
	p = list;
	while((p->list != NULL) && (p->code != code)) p++;

	/* Execute specific options */
	if(p->list){
		sprintf(msg, "vgeo_read_options(): code=%-d", code);
		check_error(execute_options(fp, p->list, ECHO) <= 0, msg);
	}
	
	/* Assign return value */
	return((char *) &code);

}

/* ------------------------------------------------------------------
	 DESCRIPTION:
	 
	    Read the grid points and allocate memory for keeping them. The value of
	    <code> specifies the method used for defining the grid points. Possible
			methods are:
			
			   GRID_LIST_METHOD: the first parameter must give the number of grid
				                   points to read as 'number=N'. There must be then N
													 successive grid points.
				 
	 RETURNS: NULL
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
char *vgeo_read_grid(FILE *fp, char *arg, int code)
#else
char *vgeo_read_grid(fp, arg, code)
FILE *fp;
char *arg;
int code;
#endif
{

	int i;

	/* Select the method for reading grid points */
	switch(code)
		{
		case GRID_LIST_METHOD:
			/* Get the number of grid points */
			check_error(fscanf(fp, " number= %d", &(opts.n_grid)) != 1, 
									"Reading number of grid points");
			/* Allocate memory */
			check_error( (opts.grid = 
										calloc((size_t) opts.n_grid, sizeof(FLOAT))) == NULL,
									 "calloc(): Grid points");
			/* Read the grids */
			for(i=0; i<(opts.n_grid); i++)
				check_error(fscanf(fp, " %f", &(opts.grid[i])) != 1,
										"Reading list of grid points");
			/* Message */
			printf("\n%% grid: method=list number=%-d:", opts.n_grid);
			for(i=0; i<(opts.n_grid); i++)
				printf("     %-8.2f\n%% ", opts.grid[i]);
			break;
		default:
			check_error(UNKNOWN_CODE, "Reading output grid");
		}

	return((char *) NULL);

}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    This function initializes pointers in specific data structures. This is
			necessary because most functions release memory from previous calls
			automatically. Therefore, some pointers have to be initialized to NULL
			before calling these the first time.

			On error, the programme will abort with a fatal exit. This is nessecary
			because most of the functions try to release allocated memory on errors.
			This would result in calling this function again.

	 RETURNS: zero
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int init_data( /* Pointer to the structure to be initialized */
							char *data, 
							/* Code specifying the structure type */
							int code, 
							/* Number of consecutive structures to initialize */
							int n
							)
#else
int ()
#endif
{

	VGEO_EXT_PRF_TYPE *vep;
	VGEO_EXT_TYPE     *ve;
	VGEO_PRFT_TYPE    *vp;
	VGEO_CTD_PRF_TYPE *vcp;
	VGEO_CTD_TYPE     *vc;
	VGEO_VGEO_TYPE    *vv;
	VGEO_TYPE         *v;
	VGEO_OPTIONS_TYPE *vo;
	VGEO_VGRID_TYPE   *vgr;

	int i, ierr=0;

	/* Check if the structure is defined */
	if(!data)
		return(0);
	
	/* Initialize structure pointers */
	for(i=0; i<n && !ierr; i++){
		switch(code)
			{
			case VGEO_EXT_PRF_CODE:
				vep		 = ((VGEO_EXT_PRF_TYPE *) data) + i;
				vep->u = vep->v = vep->z = NULL;
				vep->n = 0;
				break;
			case VGEO_EXT_CODE:
				ve = ((VGEO_EXT_TYPE *) data) + i;
				ve->data = NULL;
				ve->n = 0;
				break;
			case VGEO_PRFT_CODE:
				vp = ((VGEO_PRFT_TYPE *) data) + i;
				vp->times = NULL;
				vp->n = 0;
				break;
			case VGEO_CTD_PRF_CODE:
				vcp = ((VGEO_CTD_PRF_TYPE *) data) + i;
				vcp->rho = vcp->z   = NULL;
				vcp->n   = 0;
				break;
			case VGEO_CTD_CODE:
				vc = ((VGEO_CTD_TYPE *) data) + i;
				init_data((char *) &(vc->prev), VGEO_CTD_PRF_CODE, 1);
				init_data((char *) &(vc->this), VGEO_CTD_PRF_CODE, 1);
				break;
			case VGEO_VGEO_CODE:
				vv = ((VGEO_VGEO_TYPE *) data + i);
				vv->v = vv->z = NULL;
				vv->n = 0;
				break;
			case VGEO_CODE:
				v = ((VGEO_TYPE *) data + i);
				init_data((char *) &(v->vgeo), VGEO_VGEO_CODE,  1);
				init_data((char *) &(v->ctd) , VGEO_CTD_CODE ,  1);  
				init_data((char *) &(v->prft), VGEO_PRFT_CODE,  1);
				init_data((char *) &(v->grid), VGEO_VGRID_CODE, 1);
				break;
			case VGEO_OPTIONS_CODE:
				vo = ((VGEO_OPTIONS_TYPE *) data) + i;
				vo->grid   = NULL;
				vo->n_grid = 0;
				break;
			case VGEO_VGRID_CODE:
				vgr = ((VGEO_VGRID_TYPE *) data) + i;
				vgr->z = vgr->v = NULL;
				vgr->n = 0;
				break;
			default:
				ierr = error_found(UNDEFINED_STRUCTURE, "undefined structure\n");
				break;
			}

		/* Check for errors */
		if(ierr) exit(-1);

	}
	return(0);
}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    This function releases memory allocated for elements in specific data
			structures and initializes their pointers to NULL. The number of
			available elements is set to zero.

			On error, the programme will abort with a fatal exit. This is nessecary
			because most of the functions try to release allocated memory on errors.
			This would result in calling this function again.

	 RETURNS: zero
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int free_data(/* Pointer to data structure */
							char *data,
							/* Code for identifying the structure */
							int code,
							/* Number of consecutive structures to be deleted */
							int n
									)
#else
int ()
#endif
{

	VGEO_EXT_PRF_TYPE *vep;
	VGEO_EXT_TYPE     *ve;
	VGEO_PRFT_TYPE    *vp;
	VGEO_CTD_PRF_TYPE *vcp;
	VGEO_CTD_TYPE     *vc;
	VGEO_VGEO_TYPE    *vv;
	VGEO_TYPE         *v;
	VGEO_OPTIONS_TYPE *vo;
	VGEO_VGRID_TYPE   *vgr;

	int i, ierr=0;

	/* Check that the structure exists */
	if(!data)
		return(0);

	/* Release memory according to the code */
	for(i=0; i<n && !ierr; i++){
		switch(code)
			{
			case VGEO_EXT_PRF_CODE:
				vep = ((VGEO_EXT_PRF_TYPE *) data) + i;
				if(vep->u) free(vep->u);
				if(vep->v) free(vep->v);
				if(vep->z) free(vep->z);
				break;
			case VGEO_EXT_CODE:
				ve = ((VGEO_EXT_TYPE *) data) + i;
				free_data((char *) ve->data, VGEO_EXT_PRF_CODE, ve->n);
				if(ve->data) free(ve->data);
				break;
			case VGEO_PRFT_CODE:
				vp = ((VGEO_PRFT_TYPE *) data) + i;
				if(vp->times) free(vp->times);
				break;
			case VGEO_CTD_PRF_CODE:
				vcp = ((VGEO_CTD_PRF_TYPE *) data) + i;
				if(vcp->rho) free(vcp->rho);
				if(vcp->z)   free(vcp->z);
				break;
			case VGEO_CTD_CODE:
				vc = ((VGEO_CTD_TYPE *) data) + i;
				free_data((char *) &(vc->prev), VGEO_CTD_PRF_CODE, 1);
				free_data((char *) &(vc->this), VGEO_CTD_PRF_CODE, 1);
				break;
			case VGEO_VGEO_CODE:
				vv = ((VGEO_VGEO_TYPE *) data + i);
				if(vv->v) free(vv->v);
				if(vv->z) free(vv->z);
				break;
			case VGEO_CODE:
				v = ((VGEO_TYPE *) data + i);
				free_data((char *) &(v->vgeo), VGEO_VGEO_CODE,  1);
				free_data((char *) &(v->ctd) , VGEO_CTD_CODE ,  1);  
				free_data((char *) &(v->prft), VGEO_PRFT_CODE,  1);
				free_data((char *) &(v->grid), VGEO_VGRID_CODE, 1);
				break;
			case VGEO_OPTIONS_CODE:
				vo = ((VGEO_OPTIONS_TYPE *) data) + i;
				if(vo->grid) free(vo->grid);
				break;
			case VGEO_VGRID_CODE:
				vgr = ((VGEO_VGRID_TYPE *) data) + i;
				if(vgr->z) free(vgr->z);
				if(vgr->v) free(vgr->v);
				break;
			default:
				ierr = UNDEFINED_STRUCTURE;
				DBERROR(&ierr, "undefined structure\n");
				break;
			}
	}
	
	/* Check for errors */
	if(ierr) exit(-1);
	
	/* Initialize the structure elements */
	init_data(data, code, n);
			
	return(0);
		
}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    This function allocates memory for elements in specific data structures
			and initializes the number of available elements.

			The structure is sent to free_data() in order to release memory from
			previous calls. Hence, pointers should be initialized to NULL before
			calling this function for the first time.

			On error, already allocated memory is released and pointers are reset to
			NULL.

	 RETURNS: Error status.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int allocate_data(/* Pointer to structure to be allocated */
									char *data, 
									/* Code specifying the structure type */
									int code, 
									/* Number of elements to allocate */
									int n
									)
#else
int ()
#endif
{

	VGEO_EXT_PRF_TYPE *vep;
	VGEO_EXT_TYPE     *ve;
	VGEO_PRFT_TYPE    *vp;
	VGEO_CTD_PRF_TYPE *vcp;
	VGEO_VGEO_TYPE    *vg;
	VGEO_VGRID_TYPE   *vgr;
	
	int i, ierr=0;
	char msg[200];

	/* Check that the structure is defined */
	if(!data)
		return(0);
	
	/* Release memory from previous call */
	free_data(data, code, 1);

	/* Check that number of elements is defined */
	if(n <= 0)
		return(0);

	/* Allocate elements of structure */
	switch(code)
		{
		case VGEO_EXT_PRF_CODE:
			strcpy(msg, "External profile data structure");
			vep = (VGEO_EXT_PRF_TYPE *) data;
			vep->u = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			vep->v = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			vep->z = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			vep->n = n;
			if(!(vep->u) || !(vep->v) || !(vep->z))
				ierr = INSUFFICIENT_MEMORY;
			DBERROR(&ierr, "calloc(): External profile data");
			for(i=0; i<n && !ierr; i++)
				vep->u[i] = vep->v[i] = vep->z[i] = BADDOUBLE;
			break;
		case VGEO_EXT_CODE:
			strcpy(msg, "External data structure");
			ve = (VGEO_EXT_TYPE *) data;
			ve->data = 
				(VGEO_EXT_PRF_TYPE *) calloc((size_t) n, sizeof(VGEO_EXT_PRF_TYPE));
			ve->n = n;
			if(!(ve->data))
				ierr = INSUFFICIENT_MEMORY;
			DBERROR(&ierr, "calloc(): External data");
			init_data((char *) ve->data, VGEO_EXT_PRF_CODE, n);
			break;
		case VGEO_PRFT_CODE:
			strcpy(msg, "CTD profile times structure");
			vp = (VGEO_PRFT_TYPE *) data;
			vp->times = 
				(YMDHMS_TIME_TYPE *) calloc((size_t) n, sizeof(YMDHMS_TIME_TYPE));
			vp->n = n;
			if(!(vp->times))
				ierr = INSUFFICIENT_MEMORY;
			DBERROR(&ierr, "calloc(): VGEO profile data");
			break;
		case VGEO_CTD_PRF_CODE:
			vcp = (VGEO_CTD_PRF_TYPE *) data;
			vcp->rho = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			vcp->z   = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			if(!(vcp->rho) || !(vcp->z))
				ierr = INSUFFICIENT_MEMORY;
			DBERROR(&ierr, "calloc(): CTD profile data");
			for(i=0; i<n && !ierr; i++)
				vcp->rho[i] = vcp->z[i] = BADDOUBLE;
			vcp->n = n;
			break;
		case VGEO_VGEO_CODE:
			vg = (VGEO_VGEO_TYPE *) data;
			vg->v = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			vg->z = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			ierr = (!(vg->z) || !(vg->v)) ? INSUFFICIENT_MEMORY : 0;
			DBERROR(&ierr, "calloc(): VGEO data");
			for(i=0; i<n && !ierr; i++)
				vg->v[i] = vg->z[i] = BADDOUBLE;
			vg->n = n;
			break;
		case VGEO_VGRID_CODE:
			vgr = (VGEO_VGRID_TYPE *) data;
			vgr->z = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			vgr->v = (DOUBLE *) calloc((size_t) n, sizeof(DOUBLE));
			ierr = (!(vgr->z) || !(vgr->v)) ? INSUFFICIENT_MEMORY : 0;
			DBERROR(&ierr, "calloc(): GRID data");
			for(i=0; i<n; i++)
				vgr->z[i] = vgr->v[i] = BADDOUBLE;
			vgr->n = n;
			break;
		default:
			ierr = UNDEFINED_STRUCTURE;
			DBERROR(&ierr, "undefined structure"); 
			break;
		}
	
	/* Check for errors */
	if(ierr)
		free_data(data, code, 1);
	
	return(ierr);

}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Copy data from a structure to another one. Memory allocated for the
	    destination structure is released and reallocated to the size of the
	    source one. Hence, this structure should be initialized before calling
	    this function for the first time.

	 RETURNS: CODAS error code.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int copy_all_data(char *dst, char *src, int code)
#else
int ()
#endif
{

	VGEO_EXT_PRF_TYPE *vep0, *vep1;
	VGEO_EXT_TYPE     *ve0, *ve1;
	VGEO_CTD_PRF_TYPE *vcp0, *vcp1;

	int i, ierr=0;

	/* Initialize destination structure */
	free_data(dst, code, 1);

	switch(code)
		{
		case VGEO_EXT_PRF_CODE:
			vep0 = (VGEO_EXT_PRF_TYPE *) src;
			vep1 = (VGEO_EXT_PRF_TYPE *) dst;
			vep1->pos.lon = vep0->pos.lon;
			vep1->pos.lat = vep0->pos.lat;
			vep1->yd = vep0->yd;
			ierr = allocate_data((char *) vep1, code, vep0->n);
			DBERROR(&ierr, "allocate_data(): VGEO_EXT_PRF_TYPE");			
			for(i=0; i<(vep0->n) && !ierr; i++){
				vep1->u[i] = vep0->u[i];
				vep1->v[i] = vep0->v[i];
				vep1->z[i] = vep0->z[i];
			}
			break;
		case VGEO_EXT_CODE:
			ve0 = (VGEO_EXT_TYPE *) src;
			ve1 = (VGEO_EXT_TYPE *) dst;
			ierr = allocate_data((char *) ve1, code, ve0->n);
			DBERROR(&ierr, "allocate_data(): VGEO_EXT_TYPE");
			for(i=0; i<(ve0->n) && !ierr; i++){
				vep0 = ((VGEO_EXT_PRF_TYPE *) ve0->data) + i;
				vep1 = ((VGEO_EXT_PRF_TYPE *) ve1->data) + i;
				ierr = copy_all_data((char *) vep1, (char *) vep0, VGEO_EXT_PRF_CODE);
				DBERROR(&ierr, "copy_all_data(): VGEO_EXT_PRF_TYPE");
			}
			break;
		case VGEO_CTD_PRF_CODE:
			vcp0 = (VGEO_CTD_PRF_TYPE *) src;
			vcp1 = (VGEO_CTD_PRF_TYPE *) dst;
			vcp1->yd = vcp0->yd;
			vcp1->pos.lon = vcp0->pos.lon;
			vcp1->pos.lat = vcp0->pos.lat;
			ierr = allocate_data((char *) vcp1, VGEO_CTD_PRF_CODE, vcp0->n);
			DBERROR(&ierr, "allocate_data(): VGEO_CTD_PRF_TYPE");
			for(i=0; i<(vcp0->n) && !ierr; i++){
				vcp1->rho[i] = vcp0->rho[i];
				vcp1->z[i] = vcp0->z[i];
			}
			break;
		default:
			ierr = UNDEFINED_STRUCTURE;
			DBERROR(&ierr, "undefined structure");
			break;
		}

	return(ierr);
	
}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Used to initialize the VGEO_OPTION_TYPE and VGEO_TYPE structures before
			and after parsing of the control file. If called with <NULL>, just
			initialize default options. Else, initialize structures with current
			options.

	 RETURNS: zero
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int vgeo_init(VGEO_TYPE *vgeo)
#else
int ()
#endif
{

	FILE_NAME_TYPE new_file;
	int dba = READ_ONLY, dbm = DIR_ON_DISK;

	int i;
	char msg[200];

	/* Initialize structures with specified options */
	if(vgeo){
		/* Check consistency of reference method options */
		switch(opts.ref_method)
		{
		case ADCP_MIN_SHEAR_METHOD:
			check_error(strlen(opts.data_file) == 0, 
									"Need a file name with external data for reference");
			break;
		case CONST_DEPTH_METHOD:
			check_error(opts.z0 < 0.0, "Need a reference level");
			break;
		default:
			/* Nothing to do */
			break;
		}
		/* Check consistency of method for number of profiles */
		switch(opts.nprf_method)
			{
			case NPRF_CONST_METHOD:
				check_error(opts.n_profiles < 0,
										"Need to know the number of data profiles");
				break;
			default:
				/* Nothing to do */
				break;
			}

		/* Initialze CTD database */
		check_dbopen(CTD_DBID_CODE, ctd_db_name, dba, dbm);

		/* Initialize input files */
		vgeo->fps.n = 5;
		vgeo->fps.first = &(vgeo->fps.fp_con);
		if(strlen(opts.data_file)){
			switch(opts.data_format)
				{
				case CODAS_DB_METHOD:
					check_dbopen(ADCP_DBID_CODE, opts.data_file, dba, dbm);
					break;
				case ASCII_DATA_METHOD:
					vgeo->fps.fp_asc = check_fopen(opts.data_file, "r");
					break;
				default:
					sprintf(msg, "Unknown format for data file %s", opts.data_file);
					check_error(UNKNOWN_CODE, msg);
					break;
				}
		}

		/* Set all file pointers to zero */
		vgeo->fps.fp_con = vgeo->fps.fp_log = vgeo->fps.fp_sta =
			vgeo->fps.fp_grd = vgeo->fps.fp_asc = NULL;

		/* Contour file */
		if(opts.bits & CON_OUT){
			new_extension((char *) new_file, output, ".con");
			vgeo->fps.fp_con = check_fopen(new_file, "w");
		}
		/* Satistics file */
		if(opts.bits & STA_OUT){
			new_extension((char *) new_file, output, ".sta");
			vgeo->fps.fp_sta = check_fopen(new_file, "w");
		}
		/* Log file */
		if(opts.bits & LOG_OUT){
			new_extension((char *) new_file, output, ".log");
			vgeo->fps.fp_log = check_fopen(new_file, "w");
			set_msg_file(vgeo->fps.fp_log);
		}
		/* grid data file */
		if(opts.grid){
			new_extension((char *) new_file, (char *) output, ".grd");
			vgeo->fps.fp_grd = check_fopen(new_file, "w"); 
		}

		/* Initialize data pointers */
		init_data((char *) vgeo, VGEO_CODE, 1);

		/* Initialize regridding structure */
		if(opts.grid){
			opts.bits |= (ULONG) (to_bit(REGRID_OUTPUT_CODE));
			check_error(allocate_data((char *) &(vgeo->grid),VGEO_VGRID_CODE, 
																opts.n_grid),
									"allocate_data(): Regrid output");
			for(i=0; i<(opts.n_grid); i++)
				vgeo->grid.z[i] = opts.grid[i];
		}

		return(0);
	}
	
	/* Set default values to options */
	init_data((char *) &opts, VGEO_OPTIONS_CODE, 1);
	opts.data_format = CODAS_DB_METHOD;
	strcpy(opts.data_file, "");
	opts.data_range.type = TIME_RANGE;
	ymdhms_to_min(&(opts.data_range.ru.time.start));
	ymdhms_to_max(&(opts.data_range.ru.time.end));
	opts.bits = 0;
	opts.ref_method = CTD_MAX_DEPTH_METHOD;
	opts.z0 = BADDOUBLE;
	opts.v0 = 0.0;
	opts.data_good_bins = 15;
  opts.good_prof_ratio = 0.75;
	opts.data_interp_axis = VGEO_TIME_CODE;
	opts.nprf_method = NPRF_BCTD_METHOD;
	opts.n_profiles = -1;
	opts.adcp_uv_ref = ADCP_SHIPREF_METHOD;
	opts.ctd_sort_axis = VGEO_TIME_CODE;
	opts.ctd_steps = 1;
  opts.ctd_good_bins = 20;
	opts.ctd_time_offset = 0;
	opts.time_tolerance = 30;
	opts.vgeo_axis_dir = NORTH_DIRECTION_CODE;
	opts.vgeo_scale = 1.0;

	return(0);

}

/* ------------------------------------------------------------------
	 DESCRIPTION: 

	    Closes all databases and release memory.

	 RETURNS: CODAS error code.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int vgeo_exit(VGEO_TYPE *vgeo)
#else
int ()
#endif
{

	FILE **fp;

	int i, ierr=0;

	/* Close databases */
	ierr = dbset_cnf(CTD_DBID_CODE);
	DBERROR(&ierr, "dbset_cnf(): CTD");
	if(db->db_is_open){
		ierr = close_db();
		DBERROR(&ierr, "close_db(): CTD");
	}
	if(opts.data_format == CODAS_DB_METHOD){
		ierr = dbset_cnf(CTD_DBID_CODE);
		DBERROR(&ierr, "dbset_cnf(): ADCP");
		if(db->db_is_open){
			ierr = close_db();
			DBERROR(&ierr, "close_db(): ADCP");
		}
	}
		
	/* Close files */
	fp = vgeo->fps.first;
	for(i=0; i<(vgeo->fps.n); i++){
		if(*fp) fclose(*fp);
		fp++;
	}

	/* Release memory */
	free_data((char *) vgeo, VGEO_CODE, 1);
	free_data((char *) &opts, VGEO_OPTIONS_CODE, 1);
	
	return(ierr);

}


/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Read all CTD profiles within given time range and extract time,
	    longitude and latitude values. Then sort these profiles along the
	    specified sort axis. Stores times of the sorted profiles in the <vgeo>
	    structure.

			Memory space for the array of profile times is allocated within the
			function. If non-NULL, a call to free is performed in order to release
			memory from a previous call (typically when several time ranges are
			specified). Hence, this pointer should be initialized to NULL before
			calling this function the first time.    

	 RETURNS: Error status, or EOF if no profile founds within time range
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int init_ctd_profiles(/* Current CTD time range to consider */
											TIME_RANGE_TYPE *trng,
											/* Pointer to structure of sorted CTD profile times */
											VGEO_PRFT_TYPE *prft, 
											/* Year base to use for converting profile times into
												 decimal days. This is necessary for sorting times as
												 double data */
											int year_base    
											)
#else
int ()
#endif
{

	DOUBLE_LL_TYPE this_pos;
	YMDHMS_TIME_TYPE  this_time, *times=NULL;
	DOUBLE *lon=NULL, *lat=NULL, *yd=NULL, *ptr=NULL;

	int i, n=0, err, ierr=0;
	int *indx=NULL;

	/* Release memory from previous calls */
	free_data((char *) prft, VGEO_PRFT_CODE, 1);

	/* Select CTD database */
	ierr = dbset_cnf(CTD_DBID_CODE);
	if(DBERROR(&ierr, "dbset_cnf(): CTD\n"))
		return(ierr);

	/* Postionning to start of time range */
	ierr = search_db(TIME_SEARCH, (char *) &(trng->start));
	if(error_found(ierr, "search_db(): CTD start time\n"))
		return(ierr);
	
	/* Get the number of profiles within given time range */
	n = 0;
	while(!ierr && check_time(&this_time, &(trng->start), &(trng->end))){
		n++;
	  DBMOVE(&(opts.ctd_steps), &ierr);
	}
	/* Check for errors */
	err = db_not_eof(ierr);
	if(DBERROR(&err, "DBMOVE(): Counting CTD profiles"))
		return(err);
	if(!err && !n ) 
		return(EOF);
		
	/* Allocate space for profile index and times */
	ierr = allocate_data((char *) prft, VGEO_PRFT_CODE, n);
	if(DBERROR(&ierr, "allocate_data()"))
		return(ierr);
	indx  = (int *)              calloc((size_t) n, sizeof(int));
	times = (YMDHMS_TIME_TYPE *) calloc((size_t) n, sizeof(YMDHMS_TIME_TYPE));
	yd    = (DOUBLE *)           calloc((size_t) n, sizeof(DOUBLE));
	lon   = (DOUBLE *)           calloc((size_t) n, sizeof(DOUBLE));
	lat   = (DOUBLE *)           calloc((size_t) n, sizeof(DOUBLE));
	ierr = ( !indx || !times || !yd || !lon || !lat ) ?
		INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): CTD time and position data\n"))
		goto on_error;

	/* Move backwards to beginning of time range */
	ierr = search_db(TIME_SEARCH, (char *) &(trng->start));
	if(DBERROR(&ierr, "search_db(): Back to CTD start time\n"))
		goto on_error;

	/* Get all times and positions */
	i = 0;
	while(!ierr && (i < n)){
		/* Get time */
		ierr = get_db_time_d((DOUBLE *) yd+i, year_base);
		DBERROR(&ierr, "get_db_time_d()");
		if(!ierr){
			yd_to_ymdhmh_time(yd[i], year_base, &this_time);
			copy_ymdhms_time(times+i, &this_time);
			/* Get position */
			ierr = get_db_pos_d((DOUBLE *) &this_pos, LON_180);
			DBERROR(&ierr, "get_db_pos_d()");
		}
		if(!ierr){
			lon[i] = lon_to_lon(this_pos.lon);
			lat[i] = this_pos.lat;
			/* Move to next profile */
			DBMOVE(&(opts.ctd_steps), &ierr);
			err = db_not_eof(ierr);
			DBERROR(&err, "DBMOVE()");
		}
		if(!ierr){
			/* Update profile number */
			i++;
		}
	}
	/* Check for errors */
	ierr = db_not_eof(ierr);
	if(DBERROR(&ierr, "Getting CTD profile times and positions"))
		goto on_error;
	ierr = (i != (n-1)) ? INVALID_NUMBER : 0;
	if(DBERROR(&ierr, "Invalid profile number"))
		goto on_error;

	/* Place database at beginning of time range */
	ierr = search_db(TIME_SEARCH, (char *) &(trng->start));
	if(DBERROR(&ierr, "search_db(): Back to CTD start time"))
		goto on_error;
		 
	/* Select the sort axis */
	ptr = (DOUBLE *) NULL;
	switch(opts.ctd_sort_axis)
		{
		case VGEO_TIME_CODE:
			ptr = yd;
		case VGEO_LONGITUDE_CODE:
			ptr = lon;
			break;
		case VGEO_LATITUDE_CODE:
			ptr = lat;
			break;
		default:
			ierr = INVALID_NUMBER;
			break;
		}
	if(DBERROR(&ierr, "CTD sort axis")) 
		goto on_error;
	ierr = sort_double(n, indx, ptr);
	if(DBERROR(&ierr, "sort_double()"))
		goto on_error;
	
	/* Order times */
	for(i=0; i<n; i++)
		copy_ymdhms_time((prft->times)+i, times+indx[i]);
	
on_error:
	
	/* Release memory */
	if(indx)  free(indx);
	if(times) free(times);
	if(lon)   free(lon);
	if(lat)   free(lat);
	if(yd)    free(yd);

	if(ierr)
		free_data((char *) prft, VGEO_PRFT_CODE, 1);
		
	return(ierr);
}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Get CTD time, position, temperature and salinity and calculates density
	    and depth data for current profile. Data is regridded on the depth grid
	    specified in the P variable (e.g. 1 dbar == 1 m). Also check if the
	    number of good bins is less than what is specified in the control file.

			On error, memory allocated for CTD data is released, the corresponding
			pointers set to NULL and the number of bins is set to zero.

	 RETURNS: CODAS error status, or EOF if the number of bad data exceeds the
	          minimum allowed.

	 BUGS: Depth is calculated according to the current latitude, from the CTD
	       database, even if the use_adcp_pos flag has been set. If this value
				 is wrong, depth might not be accurate, however it is not expected to
				 be so important at present.

		     It is assumed that CTD pressure is given in decibars. To calculate
		     geostrophic current with CTD data at the same depth, regridding of
		     data is carried out to the depth values given by the pressure (e.g.
		     in the new grid, depth has the same value as the pressure, and is
		     given in meters). If CTD pressure is not given in decibars, this will
		     results in some problems.

				 It is also assumed that CTD pressure is a block variable and does not
				 vary from one profile to another. This way, all profiles will be
				 regridded on the same grid and should contain the same amount of
				 depth elements.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_ctd_data(/* CTD structure to hold density and depth data. Pointers to
										these two variables are allocated within the function. If
										non-NULL on entry, a call to free() is performed in roder
										to release memory from previous calls. Hence, they should
										be initialized to NULL before entering the function for
										the first time. */
								 VGEO_CTD_PRF_TYPE *ctd 
								 )  
#else
int ()
#endif
{

	VGEO_CTD_PRF_TYPE tmp;
	DOUBLE *sal=NULL, *tem=NULL, *p=NULL;
	DOUBLE *data=NULL;
	int i, n, nmin, ierr=0;
	double svan;

	/* Release memory from previous calls */
	free_data((char *) ctd, VGEO_CTD_PRF_CODE, 1);
	init_data((char *) &tmp, VGEO_CTD_PRF_CODE, 1);

	/* Select CTD database */
	ierr = dbset_cnf(CTD_DBID_CODE);
	if(DBERROR(&ierr, "dbset_cnf(): CTD\n"))
		return(ierr);

	/* Get the profile time and position */
	ierr = get_db_time_d(&(ctd->yd), year_base);
	if(DBERROR(&ierr, "get_db_time_d()\n"))
		return(ierr);
	ierr = get_db_pos_d( (DOUBLE *) &(ctd->pos), LON_180);
	if(DBERROR(&ierr, "get_db_pos_d\n"))
		return(ierr);
	
	/* Get CTD data */
	nmin = BADINT;
	ierr = get_db_data_d("TEMPERATURE", "", &tem, (unsigned int *) &n, year_base);
	if(DBERROR(&ierr, "get_db_data_d(): TEMPERATURE\n"))
		goto on_error;
	nmin = min_val(nmin, n);
	ierr = get_db_data_d("SALINITY", "", &sal, (unsigned int *) &n, year_base);
	if(DBERROR(&ierr, "get_db_data_d(): SALINITY\n"))
		goto on_error;
	nmin = min_val(nmin, n);
	ierr = get_db_data_d("P", "", &p, (unsigned int *) &n, year_base);
	if(DBERROR(&ierr, "get_db_data_d(): P\n"))
		goto on_error;
	nmin = min_val(nmin, n);
	
	/* Check the number of data */
	ierr = ( (nmin == BADINT) || (nmin <= 0) ) ? INVALID_NUMBER : 0;
	if(ierr) goto on_error;
	
	/* Allocate memory for density and actual depth */
	ierr = allocate_data((char *) ctd, VGEO_CTD_PRF_CODE, nmin);
	if(DBERROR(&ierr, "allocate_data(): CTD density and depth\n"))
		goto on_error;
	ierr = allocate_data((char *) &tmp, VGEO_CTD_PRF_CODE, nmin);
	if(DBERROR(&ierr, "allocate_data(): TMP density and depth\n"))
		goto on_error;

	/* Check that time and positions are defined */
	ierr = (ctd->yd > ADJ_BADFLOAT) ? INVALID_TIME : 0;
	ierr = ((ctd->pos.lon > ADJ_BADFLOAT) 
					|| (ctd->pos.lat > ADJ_BADFLOAT)) ? INVALID_POSITION : 0;
	if(ierr) goto on_error;
 
	/* Calculate actual depth and density */
	for(i=0; i<(ctd->n); i++){
		tmp.rho[i] = BADDOUBLE;
		tmp.z[i] = (DOUBLE) p_depth(p[i], tmp.pos.lat);
		ctd->z[i] = p[i];
		if( (sal[i] < ADJ_BADFLOAT) && (tem[i] < ADJ_BADFLOAT) ){
			svan = SVAN(sal[i], tem[i], p[i], tmp.rho+i);
			tmp.rho[i] += 1000.0;
		}
	}

	/* Regrid density on pressure data, filling gaps */
	n = s_regrid_d(tmp.z, tmp.rho, ctd->z, ctd->rho, tmp.n, ctd->n);
	ierr = (n == EOF) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "s_regrid_d(): regrid CTD"))
		goto on_error;

	/* Get the number of good bins */
	n = 0;
	for(i=0; i<(ctd->n); i++){
		if(ctd->rho[i] < ADJ_BADFLOAT)
			n++;
	}
	ierr = (n < opts.ctd_good_bins) ? EOF : 0;
	if(ierr) goto on_error;

on_error:

	/* Release memory */
	if(data) free(data);
	if(tem)  free(tem);
	if(sal)  free(sal);
	if(p)    free(p);
	free_data((char *) &tmp, VGEO_CTD_PRF_CODE, 1);

	/* Check error */
	if(ierr)
		free_data((char *) ctd, VGEO_CTD_PRF_CODE, 1);

	return(ierr);

}

/* ------------------------------------------------------------------
	 DESCRIPTION:
	 
	    Move to next CTD profile. If <prft> is defined, moves to the profile
	    whose time is the first element of this list. If undefined, moves
	    <ctd_steps> profiles forward, as specified in the option structure.

	    On output, the list of profile times shifted so that its first element
	    is the time for the next profile. The number of entries in this list,
	    <n>, is also updated, and the list is filled with bad times at the end.
	 
			RETURNS: CODAS error code. If using a list of profile times,
			         SEARCH_BEYOND_END will be returned if no more profile times are
			         available
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int next_ctd_profile(/* Pointer to array of CTD profile times */
										 VGEO_PRFT_TYPE *prft
										 )
#else
int ()
#endif
{

	YMDHMS_TIME_TYPE this_time, *t;

	LONG tol;
	int i, err, ierr=0;

	/* Select CTD database */
	ierr = dbset_cnf(CTD_DBID_CODE);
	if(DBERROR(&ierr, "dbset_cnf(): CTD\n")) 
		return(ierr);

	if(!(prft->times)){
		/* Move the required number of steps */
		DBMOVE(&(opts.ctd_steps), &ierr);
		err = db_not_eof(ierr);
		DBERROR(&err, "DBMOVE(): CTD");
		return(ierr);
	}
	
	/* Check for available times */
	if(prft->n <= 0)
		return(SEARCH_BEYOND_END);
	ierr = (!(prft->times)) ? UNDEFINED_STRUCTURE : 0;
	if(DBERROR(&ierr, "Undefined CTD profile times"))
		return(ierr);
	
	/* Set the time for searching just before the next profile time */
	tol = (LONG) (opts.time_tolerance) * (-1);
	DIFTIM(prft->times, &this_time, (LONG *) &tol);

	/* Position to profile */
	ierr = search_db(TIME_SEARCH, (char *) &this_time);
	if(DBERROR(&ierr, "search_db()"))
		return(ierr);
		
	/* Update the list */
	for(i=1; i<(prft->n); i++){
		t = prft->times + i;
		copy_ymdhms_time(t-1, t);
	}
	t = prft->times + (prft->n) - 1;
	ymdhms_to_bad(t);
	(prft->n)--;
		
	return(ierr);
}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Get ADCP currents from a CODAS database to calculate the reference for
	    the current geostrophic profile. If requested, ADCP currents are flagged
	    using the bits in the PROFILE_FLAGS variable, and converted to absolute
	    currents, using the specified reference method.

	    If an error occurs, allocated memory is released before leaving the
	    function, pointers set to NULL, and number of profiles set to zero.

	 RETURNS: CODAS error code. If an ADCP profile cannot be found, returns
			      SEARCH_BEYOND_END.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_adcp_db_data( /* Structure to hold all profiles of external ADCP data.
												 Memory space is allocated within the function. If
												 non_NULL on entry, a call to free_ext_data() is
												 performed in order to release memory allocated from
												 previous calls. Hence, this posinter should be
												 initialized to NULL before entering this function the
												 first time. */
										 VGEO_EXT_TYPE *ext, 
										 /* CTD profile values from the two current consecutive
												profiles. */
										 VGEO_CTD_TYPE *ctd,
										 /* Geostrophic current structure */
										 VGEO_VGEO_TYPE *vg,
										 /* Year base for time conversion into decimal days */
										 int year_base
										 )
#else
int ()
#endif
{

	VGEO_EXT_PRF_TYPE *d;

	UBYTE *pf=NULL;
	DOUBLE *ship_data=NULL, u_ship, v_ship, ship_dir, ship_speed;

	int i, j, err, ierr=0;
	int nmin, nd, id;
	int n0, n1, offset, nprof, steps=1;

	/* Release memory from previous calls */
	free_data((char *) ext, VGEO_EXT_CODE, 1);

	/* Select ADCP database */
	ierr = dbset_cnf(ADCP_DBID_CODE);
	if(DBERROR(&ierr, "dbset_cnf(): ADCP"))
		return(ierr);
	
	/* Get the profile offset and number of profiles to get from this database
		 position, according to the selected method */
	switch(opts.nprf_method)
		{
		case NPRF_CONST_METHOD:
			/* Count number of profiles to vgeo position */
			n0 = -1;
			switch(opts.data_interp_axis)
				{
				case VGEO_TIME_CODE:
					ierr = seek_db((char *) &(vg->yd), DAY_SEEK, 
												 &(opts.data_range), year_base, &n0);
					break;
				case VGEO_LONGITUDE_CODE:
					ierr = seek_db((char *) &(vg->pos.lon), LONGITUDE_SEEK, 
											&(opts.data_range), LON_180, &n0);
					break;
				case VGEO_LATITUDE_CODE:
					ierr = seek_db((char *) &(vg->pos.lat), LATITUDE_SEEK, 
											 &(opts.data_range), 0, &n0);
					break;
				default:
					ierr = error_found(INVALID_NUMBER, "interpolation axis\n");
					break;
				}
			/* Check for errors */
			ierr = (ierr == NO_SUCH_BLOCK_PROFILE) ? SEARCH_BEYOND_END : ierr;
			err = db_not_eof(ierr);
			DBERROR(&err, "seek_db(): vgeo profile");
			if(ierr) return(ierr);
			/* Calculate offset and number of profiles */				
			offset = n0 - ( (int) (opts.n_profiles / 2) );
			nprof = opts.n_profiles;
			break;
		case NPRF_BCTD_METHOD:
			/* Move to first CTD profile */
			n0 = -1;
			switch(opts.data_interp_axis)
				{
				case VGEO_TIME_CODE:
					ierr = seek_db((char *) &(ctd->prev.yd), DAY_SEEK, 
												 &(opts.data_range), year_base, &n0);
					break;
				case VGEO_LONGITUDE_CODE:
					ierr = seek_db((char *) &(ctd->prev.pos.lon), LONGITUDE_SEEK, 
												 &(opts.data_range), LON_180, &n0);
					break;
				case VGEO_LATITUDE_CODE:
					ierr = seek_db((char *) &(ctd->prev.pos.lat), LATITUDE_SEEK, 
												 &(opts.data_range), 0, &n0);
					break;
				default:
					ierr = error_found(INVALID_NUMBER, "interpolation axis\n");
					break;
				}
			/* Check for errors */
			ierr = (ierr == NO_SUCH_BLOCK_PROFILE) ? SEARCH_BEYOND_END : ierr;
			err = db_not_eof(ierr);
			DBERROR(&err, "seek_db(): First CTD");
			if(ierr) return(ierr);
			/* Move to second CTD profile */
			n1 = -1;
			switch(opts.data_interp_axis)
				{
				case VGEO_TIME_CODE:
					ierr = seek_db((char *) &(ctd->this.yd), DAY_SEEK, 
												 &(opts.data_range), year_base, &n1);
					break;
				case VGEO_LONGITUDE_CODE:
					ierr = seek_db((char *) &(ctd->this.pos.lon), LONGITUDE_SEEK, 
												 &(opts.data_range), LON_180, &n1);
					break;
				case VGEO_LATITUDE_CODE:
					ierr = seek_db((char *) &(ctd->this.pos.lat), LATITUDE_SEEK, 
												 &(opts.data_range), 0, &n1);
					break;
				default:
					ierr = error_found(INVALID_NUMBER, "interpolation axis\n");
					break;
				}
			/* Check for errors */
			ierr = (ierr == NO_SUCH_BLOCK_PROFILE) ? SEARCH_BEYOND_END : ierr;
			err = db_not_eof(ierr);
			DBERROR(&err, "seek_db(): Second CTD");
			if(ierr) return(ierr);
			/* Calculate offset and number of profiles */
			offset = min_val(n0, n1);
			nprof  = abs_val(n1-n0) + 1;
			break;
		default:
			ierr = error_found(INVALID_NUMBER, "profile method\n");
			break;
		}
	/* Check for errors */
	if(DBERROR(&ierr, "profile offset and number"))
		return(ierr);

	/* Position DB at beginning of range */
	ierr = goto_start_of_range(&(opts.data_range));
	ierr = db_not_beof(ierr);
	if(DBERROR(&ierr, "goto_start_of_range()"))
		return(SEARCH_ERROR);

	/* Move offset profiles */
	DBMOVE(&offset, &ierr);
	if(DBERROR(&ierr, "DBMOVE()"))
		return(ierr);
	
	/* Allocate space for external data */
	ierr = allocate_data((char *) ext, VGEO_EXT_CODE, nprof);
	if(DBERROR(&ierr, "allocate_data(): external data"))
		return(ierr);

	/* Get profile data */
	for(i=0; i<(ext->n) && !ierr; i++){
		d = (ext->data) + i;
		d->n = 0;
		/* Get time */
		ierr = get_db_time_d(&(d->yd), year_base);
		DBERROR(&ierr, "get_db_time_d()\n");
		/* Get position */
		if(!ierr){
			ierr = get_db_pos_d((DOUBLE *) &(d->pos), LON_180);
			DBERROR(&ierr, "get_db_pos_d()\n");
		}
		/* Get data */
		nmin = BADINT;
		if(!ierr){
			ierr = get_db_data_d("DEPTH", "", &(d->z), (unsigned int *) &nd, 0);
			DBERROR(&ierr, "get_db_data_f(): DEPTH\n");
		}
		if(!ierr){
			nmin = min_val(nmin, nd);
			ierr = get_db_data_d("U", "", &(d->u), (unsigned int *) &nd, 0);
			DBERROR(&ierr, "get_db_data_f(): U\n");
		}
		if(!ierr){
			nmin = min_val(nmin, nd);
			ierr = get_db_data_d("V", "", &(d->v), (unsigned int *) &nd, 0);
			DBERROR(&ierr, "get_db_data_f(): V\n");
		}
		if(!ierr)
			nmin = min_val(nmin, nd);
		/* Flag data using the profile flags of current profile */
		if(!ierr && FLAG_DATA){
			ierr = get_db_data("PROFILE_FLAGS", "", (char **) &pf, &id, 
												 (unsigned int *) &nd);
			DBERROR(&ierr, "get_db_data_f(): PROFILE_FLAGS\n");
			if(!ierr){
				nmin = min_val(nmin, nd);
				for(j=0; j<nmin; j++){
					if(pf[j] & ALL_BITS){
						(d->u)[j] = BADDOUBLE;
						(d->v)[j] = BADDOUBLE;
					}
				}
			}
		}
		/* Update the number of common data and calculate absolute currents */
		if(!ierr){
			d->n = nmin;
			/* Calculate absolute current */
			switch(opts.adcp_uv_ref)
				{
				case ADCP_NOREF_METHOD:
					u_ship = v_ship = 0.0;
					break;
				case ADCP_NAVREF_METHOD:
					ierr = get_db_data_d("NAVIGATION", "speed", &ship_data, 
															 (unsigned int *) &nd, 0);
					err = db_not_eof(ierr);
					DBERROR(&err, "get_db_data(): ship speed\n");
					ship_speed = *ship_data;
					ierr = get_db_data_d("NAVIGATION", "direction", &ship_data,
															 (unsigned int *) &nd, 0);
					err = db_not_eof(ierr);
					DBERROR(&err, "get_db_data(): ship direction\n");
					ship_dir = *ship_data;
					/* Calculate components of ship velocity */
					u_ship = KNOT_TO_MPS * ship_speed * cos(to_rad(90.0 - ship_dir));
					v_ship = KNOT_TO_MPS * ship_speed * sin(to_rad(90.0 - ship_dir));
					break;
				case ADCP_SHIPREF_METHOD:
					/* Get ship velocity */
					ierr = get_db_data_d("ACCESS_VARIABLES", "U_ship_absolute", 
															 &ship_data, (unsigned int *) &nd, 0);
					err = db_not_eof(ierr);
					DBERROR(&err, "get_db_data(): U_ship_absolute\n");
					u_ship = *ship_data;
					ierr = get_db_data_d("ACCESS_VARIABLES", "V_ship_absolute", 
															 &ship_data, (unsigned int *) &nd, 0);
					err = db_not_eof(ierr);
					DBERROR(&err, "get_db_data(): V_ship_absolute\n");
					v_ship = *ship_data;
					break;
				default:
					ierr = INVALID_NUMBER;
					DBERROR(&ierr, "reference method\n");
					break;
				}
			/* Loop through bins */
			for(j=0; j<(d->n) && !ierr; j++){
				if((d->u[j] < ADJ_BADFLOAT) && (u_ship < ADJ_BADFLOAT))
					d->u[j] += u_ship;
				else
					d->u[j] = BADDOUBLE;
				if((d->v[j] < ADJ_BADFLOAT) && (v_ship < ADJ_BADFLOAT))
					d->v[j] += v_ship;
				else
					d->v[j] = BADDOUBLE;
			}

		}
		/* Move to next profile */
		if(!ierr){
			DBMOVE(&steps, &ierr);
			DBERROR(&ierr, "DBMOVE()");
		}
	}
	/* Check for errors */
	if(ierr) goto on_error;

on_error:
		
	/* Release memory */
	if(pf)  free(pf);
	if(ship_data) free(ship_data);
	
	/* Check for errors */
	if(ierr)
		free_data((char *) ext, VGEO_EXT_CODE, 1);

	return(ierr);
			
}
	 
/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Get external data according to the specified file format and profile
			retrieval method. It removes in addtion bins of external data that are
			not within the geostrophic current depth range (this to assure that the
			reference level will be within the geostrophic depth range). Elements in
			bins that do contain at least one bad value are all set to bad (except
			depth which is needed to reference bins). Furthermore, profiles with
			less than the specified good bins parameter are deleted, so they won't
			be taken into account for calculating the reference. Finally, checks if
			the ratio of good profiles is satisfactory, as specified in the control
			file.

			On error, memory allocated for data is released, its pointer set to NULL
			and the number of profiles set to zero.

	 RETURNS: CODAS error code. If unable to get all necessary data, returns
	          SEARCH_BEYOND_END to indicate the end of accessible data, or that
	          a too low amount of good profiles has been read.

   BUGS: It is assumed that the DEPTH values are identical for all profiles
         (since it is a BLOCK_VARIABLE, this should be the case). It is also
         assumed that no depth values are set to bad.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_external_data(/* Pointer to array of external data. Memory space for
												 data is allocated within memory. If non_NULL on
												 entry, a call to free_ext_data() is performed in
												 order to release memory from previous calls. Hence,
												 the array should be initialized to NULL before
												 calling this function for the first time. */
											VGEO_EXT_TYPE *ext,
											/* Pointer to the two conscutive CTD data profiles used
												 to calculate current geostrophic current */
											VGEO_CTD_TYPE *ctd, 
											/* Pointer to the geostrophic profile data */
											VGEO_VGEO_TYPE *vg)
#else
int ()
#endif
{

	VGEO_EXT_TYPE     tmp;
	VGEO_EXT_PRF_TYPE *d, *d0, *d1;

	double zmin, zmax;
	int err, ierr=0;
	int i, j;
	int good_p, good_b;
	int *indx;

	/* Release memory from previous calls */
	free_data((char *) ext, VGEO_EXT_CODE, 1);
	init_data((char *) &tmp, VGEO_EXT_CODE, 1);

	/* Get external data according to the selected method */
	switch(opts.data_format)
		{
		case ASCII_DATA_METHOD:
			ierr = error_found(INVALID_NUMBER, 
												 "Sorry: Ref data from ASCII files not available\n");
			break;
		case CODAS_DB_METHOD:
			ierr = get_adcp_db_data(&tmp, ctd, vg, year_base);
			err = db_not_eof(ierr);
			DBERROR(&err, "get_adcp_db_data()\n");
			break;
		default:
			ierr = error_found(INVALID_NUMBER, "data format\n");
			break;
		}
	/* Check for errors */
	if(ierr) return(ierr);
	
	/* Find min and max depths of geostrophic current */
	zmin = BADDOUBLE;
	zmax = -BADDOUBLE;
	for(i=0; i<(vg->n); i++){
		zmin = min_val(zmin, vg->z[i]);
		zmax = max_val(zmax, vg->z[i]);
	}

	/* Allocate memory for index array to good profiles */
	indx = (int *) calloc((size_t) tmp.n, sizeof(int));
	ierr = (!indx) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc(): Index to good profiles\n"))
		goto on_error;

	/* Flag data and count good profiles */
	good_p = 0;
	for(i=0; i<(tmp.n); i++){
		d = tmp.data + i;
		good_b = 0;
		for(j=0; j<(d->n); j++){
			/* Check if within vgeo depth range */
			if((d->z[j] <= zmin) || (d->z[j] >= zmax))
				d->u[j] = d->v[j] = BADFLOAT;
			/* Cound good bins */
			if((d->u[j] < ADJ_BADFLOAT) && (d->v[j] < ADJ_BADFLOAT))
				good_b++;
			else
				d->u[j] = d->v[j] = BADDOUBLE;
		}
		/* Check the number of good bins */
		if(good_b >= opts.data_good_bins){
			indx[good_p] = i;
			good_p++;
		}
	}
	/* Check the number of good profiles */
	ierr = (good_p < ((int) (opts.good_prof_ratio * tmp.n))) ? 
		SEARCH_BEYOND_END : 0;
	err = db_not_eof(ierr);
	if(DBERROR(&err, "Good profiles\n") || ierr)
		goto on_error;

	/* Allocate memory for good profiles only */
	ierr = allocate_data((char *) ext, VGEO_EXT_CODE, good_p);
	if(DBERROR(&ierr, "allocate_data(): good profiles\n"))
		goto on_error;
	
	/* Copy good profiles */
	for(i=0; i<(ext->n) && !ierr; i++){
		d0 = tmp.data + indx[i];
		d1 = (ext->data) + i;
		ierr = allocate_data((char *) d1, VGEO_EXT_PRF_CODE, d0->n);
		DBERROR(&ierr, "allocate_data(): copy for profile\n");
		if(!ierr){
			ierr = copy_all_data((char *) d1, (char *) d0, VGEO_EXT_PRF_CODE);
			DBERROR(&ierr, "copy_all_data(): good profile\n");
		}
	}
	if(ierr) goto on_error;
	
on_error:

	/* Release memory */
	if(indx) free(indx);
	if(tmp.data)
		free_data((char *) &tmp, VGEO_EXT_CODE, 1);
	
	/* Check for errors */
	if(ierr)
		free_data((char *) ext, VGEO_EXT_CODE, 1);

	return(ierr);

}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Set the reference according to the specified method.

	 RETURNS: CODAS error code
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int get_reference(/* Pointer to global structure */
									VGEO_TYPE *vgeo
									)
#else
int ()
#endif
{

	VGEO_EXT_TYPE     ext;
	VGEO_EXT_PRF_TYPE *d, *d0, *d1;

	VGEO_REF_TYPE *ref;
	VGEO_CTD_TYPE *ctd;
	VGEO_VGEO_TYPE *vg;

	UNISTAT_TYPE s, u;

	double *u0=NULL, *u1=NULL, *z=NULL, *sh=NULL;
	double ds;
	int err, ierr=0;
	int j, i;
	int n_bins, imin;

	/* Assign pointers */
	init_data((char *) &ext, VGEO_EXT_CODE, 1);
	ref = &(vgeo->ref);
	ctd = &(vgeo->ctd);
	vg  = &(vgeo->vgeo);

	/* Initialize the reference structure */
	ref->z = ref->v = ref->e = BADDOUBLE;
	ref->n = 0;

	/* Select the reference method */
	switch(opts.ref_method)
		{
		case CONST_DEPTH_METHOD:
			/* Copy options into reference structure */
			ref->z = opts.z0;
			ref->v = opts.v0;
			ref->e = 0.0;
			ref->n = 0;
			break;
		case CTD_MAX_DEPTH_METHOD:
			/* Get the maximal VGEO good depth. This should correspond the maximal
				 common CTD depth. It is an error if this is the first bin, because we
				 should have at least several good bins within a geostrophic profile */
			i = vg->n - 1;
			while((vg->z[i] > ADJ_BADFLOAT) && (i > 0)) i--;
			if(i == 0)
				ref->z = BADDOUBLE;
			else
				ref->z = vg->z[i];
			ref->v = opts.v0;
			ref->e = 0.0;
			ref->n = 0;
			break;
		case ADCP_MIN_SHEAR_METHOD:
			/* Get external adcp data. This function must prepare consecutive
				 profiles of external data, where each profile can be used in the
				 calculations below. */
			ierr = get_external_data(&ext, ctd, vg);
			err = db_not_eof(ierr);
			DBERROR(&err, "get_external_data()\n");
			if(ierr) break;
			if(!ierr){
				/* Get the maximal number of common bins */
				n_bins = MAXINT;
				for(i=0; i<(ext.n); i++){
					d = (ext.data) + i;
					n_bins = min_val(n_bins, d->n);
				}
				/* Allocate space for current projection and depth */
				u0= (double *) calloc((size_t) n_bins, sizeof(double));
				u1= (double *) calloc((size_t) n_bins, sizeof(double));
				z = (double *) calloc((size_t) n_bins, sizeof(double));
				sh= (double *) calloc((size_t) n_bins, sizeof(double));
				ierr = (!u0 || !u1 || !z || !sh) ? INSUFFICIENT_MEMORY : 0;
				DBERROR(&ierr, "calloc(): Minimal shear data\n");
			}
			if(!ierr){
				/* Fill the depth array */
				for(i=0; i<n_bins; i++)
					z[i] = BADDOUBLE;
				for(i=0; i<(ext.n); i++){
					d = (ext.data) + i;
					for(j=0; j<n_bins; j++)
						if(d->z[j] < ADJ_BADFLOAT)
							z[j] = d->z[j];
				}
				
				/* Initialize statistics */
				allocate_unistat(&s, n_bins, "Horizontal Shear",
												 U_VARIANCE);
				allocate_unistat(&u, n_bins, "Projected Current",
												 U_VARIANCE);
				zero_unistat(&s);
				zero_unistat(&u);
				/* Initialize first profile with current projection */
				d0 = (ext.data);
				for(i=0; i<n_bins; i++){
					u0[i] = BADDOUBLE;
					if(z[i] < ADJ_BADFLOAT)
						u0[i] = (vg->unit_vector[0] * d0->u[i]) + 
							(vg->unit_vector[1] * d0->v[i]);
				}
				update_unistat_d(&u, u0, n_bins); 
				/* Loop through the next profiles */
				for(i=1; i<(ext.n); i++){
					d0 = (ext.data) + i - 1;
					d1 = (ext.data) + i;
					/* Calculate distance bewteen profiles */
					ds = posdiff_meter_d(d0->pos.lon, d0->pos.lat, 
															 d1->pos.lon, d1->pos.lat);
					/* Calculate current projection of current profile and horizontal
						 shear values */
					for(j=0; j<n_bins; j++){
						u1[j] = sh[j] = BADDOUBLE;
						if(z[j] < ADJ_BADFLOAT)
							u1[j] = (vg->unit_vector[0] * d1->u[j])
								* (vg->unit_vector[1] * d1->v[j]);
						if((ds > MINFLOAT) 
							 && (u0[j] < ADJ_BADFLOAT)
							 && (u1[j] < ADJ_BADFLOAT))
							sh[j] = (u1[j] - u0[j])/ds;
						/* Copy current projection in previous profile */
						u0[j] = u1[j];
					}
					/* Update statistics */
					update_unistat_d(&s, sh, n_bins);
					update_unistat_d(&u, u1, n_bins);
				}
				/* Calculate statistics */
				calculate_unistat(&s);
				/* Find the level with minimum shear calculated with a sufficient
					 number of points */
				imin = BADINT;
				ref->e = BADDOUBLE;
				for(i=0; i<n_bins; i++){
					if((s.sum[i] < ADJ_BADFLOAT)
						 && (z[i] < ADJ_BADFLOAT)
						 && (opts.good_prof_ratio 
								 <=  ((double) s.npts[i])/((double) (ext.n)))
						 && (abs_val(s.sum[i]) < ref->e)){
						ref->e = abs_val(s.sum[i]);
						imin = i;
					}
				}
				if(error_found( imin == BADINT, "Reference not found\n"))
					ierr = NO_SUCH_ELEMENT;
			}
			if(!ierr){
				/* Fill the reference structure */
				ref->z = z[imin];
				ref->v = u.sum[imin];
				ref->e = sqrt(u.sumsq[imin]);
				ref->n = u.npts[imin];
			}
			/* Release memory used to calculate statistics */
			free_unistat(&s);
			free_unistat(&u);
			break;
		default:
			ierr = error_found(INVALID_NUMBER, "Reference method\n");
			break;
		} 
	if(ierr) goto on_error;
	
on_error:

	/* Release memory */
	if(u0) free(u0);
	if(u1) free(u1);
	if(sh) free(sh);
	if(z) free(z);
	if(ext.data) free_data((char *) &ext, VGEO_EXT_CODE, 1);
	
	return(ierr);

}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	 Calculate the vertical derivative of the geostrophic current, and the unit
	 vector defining its direction. Both CTD profile structures <this> and
	 <prev> must be filled. Memory space for geostrophic data is allocated
	 within the function. On input, if the pointer to geostrophic data are
	 non_NULL on entry, a call to free() is performed in order to release memory
	 from a previous call. Hence, these should be initialized to NULL before
	 calling this function for the first time. The function does also checks is
	 the final number of geostrophic data is at least the number specified in
	 the control file.

	 On error, memory allocated for geostrophic data is released, the
	 corresponding pointers set to NULL, and the number of bins is set to zero.

	 RETURNS: CODAS error code, or EOF if there were not enough good bins.
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
int calculate_vgeo(/* Pointer to structure for geostrophic data. Required
											memory for these data is allocated within the function.
											If non-NULL on entry, allocated space from previous
											call is released. Hence, these should be initialized to
											NULL before calling this function the first time */
									 VGEO_VGEO_TYPE *vg, 
									 /* Pointer to two consecutive CTD profile data. It is
											assumed that these have been regridded on a common grid.
											*/ 
									 VGEO_CTD_TYPE *ctd)
#else
int ()
#endif
{

	VGEO_CTD_PRF_TYPE *c0, *c1;

	double f, ds, dx, dy;
	int ierr=0;
	int n, i;

	/* Initialize pointers */
	c0 = &(ctd->prev);
	c1 = &(ctd->this);

	/* Initializations */
	free_data((char *) vg, VGEO_VGEO_CODE, 1);

	/* Swap CTD profiles according to selected direction */
	switch(opts.vgeo_axis_dir)
		{
		case NORTH_DIRECTION_CODE:
			if(c1->pos.lat < c0->pos.lat){
				c0 = &(ctd->this);
				c1 = &(ctd->prev);
			}
			break;
		case SOUTH_DIRECTION_CODE:
			if(c1->pos.lat > c0->pos.lat){
				c0 = &(ctd->this);
				c1 = &(ctd->prev);
			}
			break;
		case WEST_DIRECTION_CODE:
			if(c1->pos.lon > c0->pos.lon){
				c0 = &(ctd->this);
				c1 = &(ctd->prev);
			}
			break;
		case EAST_DIRECTION_CODE:
			if(c1->pos.lon < c0->pos.lon){
				c0 = &(ctd->this);
				c1 = &(ctd->prev);
			}
			break;
		default:
			ierr = error_found(INVALID_NUMBER, "axis direction\n");
			break;
		}
	if(ierr) return(ierr);

	/* Allocate memory for the minimum number of of depth levels */
	ierr = allocate_data((char *) vg, VGEO_VGEO_CODE, min_val(c0->n, c1->n));
	if(DBERROR(&ierr, "allocate_data(): vg"))
		goto on_error;

	/* Calculate time and position */
	vg->pos.lat = (c0->pos.lat + c1->pos.lat)/2.0;
	vg->pos.lon = (c0->pos.lon + c1->pos.lon)/2.0;
	vg->yd      = (c0->yd + c1->yd)/2.0;
	
	/* Correct CTD time */
	vg->yd += ((DOUBLE) (opts.ctd_time_offset))/86400.0;

	/* Coriolis Parameter */
	f = 4.0 * M_PI * sin(M_PI*vg->pos.lat/180.0) / 86400.0;
	
	/* Calculate distance */
	ds = posdiff_meter_d(c0->pos.lon, c0->pos.lat, c1->pos.lon, c1->pos.lat);
	ierr = (abs_val(ds) > MINFLOAT) ? 0 : INVALID_NUMBER;
	if(DBERROR(&ierr, "zero distance between CTD"))
		goto on_error;

	/* Calculate unit vector. This should correspond to the direction othogonal
		 to the ship track, positive startboard of the chosen direction axis */
	dx = (double) dlon_to_meters(      c1->pos.lon - c0->pos.lon,
																0.5*(c1->pos.lat + c0->pos.lat));
	dy = (double) dlat_to_meters(      c1->pos.lat - c0->pos.lat,
																0.5*(c1->pos.lat + c0->pos.lat));
	vg->unit_vector[0] = dy/ds;
	vg->unit_vector[1] = -dx/ds;

	/* Loop through the depth levels */
	n = 0;
	for(i=0; i<(vg->n); i++){
		vg->v[i] = BADDOUBLE;
		vg->z[i] = c0->z[i];
		if((c0->rho[i] < ADJ_BADFLOAT) && (c1->rho[i] < ADJ_BADFLOAT)){
			/* Calculate vertical derivative of geostrophic current */
			vg->v[i] = - GRAV(vg->z[i], vg->pos.lat) /
				f*(1.0 - c0->rho[i]/c1->rho[i])/ds; 
			n++;
		}
	}

	/* Check that there are enough good bins */
	ierr = (n < opts.ctd_good_bins) ? EOF : 0;

on_error:

	/* Check for errors */
	if(ierr)
		free_data((char *) vg, VGEO_VGEO_CODE, 1);
	
	return(ierr);
	
}

/* ------------------------------------------------------------------
	 DESCRIPTION:

	    Integrate vertically geostrophic current, using the given reference. On
	    output, the geostrophic depths are updated to the middle integration
	    segments, and geostrophic current corresponds to absolute values.
			
	 RETURNS: CODAS error code or EOF if reference is undefined.
	 ------------------------------------------------------------------ */
#if PROTOTYPE_ALLOWED
int integrate_vgeo(/* Pointer to vertical derivative of geostrophic current */
									 VGEO_VGEO_TYPE *vg, 
									 /* Pointer to reference level structure */
									 VGEO_REF_TYPE *ref)
#else
int ()
#endif
{

	DOUBLE *v, *z;
	int *indx;

	DOUBLE v0;
	int ierr=0;
	int i, j, n, j0, j1;

	/* Control if reference level is defined */
	if((ref->z > ADJ_BADFLOAT) || (ref->v > ADJ_BADFLOAT))
		return(EOF);
	
	/* Allocate memory for temprary data and index to good bins */
	v = (DOUBLE *) calloc((size_t) vg->n, sizeof(DOUBLE));
	z = (DOUBLE *) calloc((size_t) vg->n, sizeof(DOUBLE));
	indx = (int *) calloc((size_t) vg->n, sizeof(int));
	ierr = (!v || !z || !indx) ? INSUFFICIENT_MEMORY : 0;
	if(DBERROR(&ierr, "calloc()"))
		goto on_error;

	/* Initialize arrays */
	j = 0;
	for(i=0; i<(vg->n); i++){
		v[i] = BADDOUBLE;
		if((vg->v[i] < ADJ_BADFLOAT) && (vg->z[i] < ADJ_BADFLOAT)){
			indx[j] = i;
			j++;
		}
	}
	n = j;
	ierr = (j <= 2) ? INVALID_NUMBER : 0;
	if(DBERROR(&ierr, "not enough data"))
		goto on_error;
	
	/* Initial value */
	v[0] = 0.0;
	z[0] = vg->z[indx[0]];
	for(i=1; i<n; i++){
		j0 = indx[i-1];
		j1 = indx[i];
		z[i] = 0.5 * (vg->z[j0] + vg->z[j1]);
		v[i] = v[i-1] + 0.5 * (vg->v[j1] + vg->v[j0]) * (vg->z[j1] - vg->z[j0]);
	}

	/* Get integral value from first bin to reference level */
	i = 1;
	while((i < n) && (vg->z[indx[i]] < ref->z)) i++;
	ierr = (i < n) ? 0 : INVALID_NUMBER;
	if(DBERROR(&ierr, "Reference depth not found!"))
		goto on_error;
	j0 = indx[i-1];
	j1 = indx[i];
	v0 = (ref->z - vg->z[j0])*(vg->v[j1] - vg->v[j0])/(vg->z[j1] - vg->z[j0]);
	v0 = v[i-1] + 0.5 * (vg->v[j0] + v0) * (ref->z - vg->z[j0]);
				
	/* Calculate integral relative to reference depth */
	for(i=0; i<n; i++)
		v[i] -= v0;
	
	/* Calculate absolute values and backup data */
	for(i=0; i<n; i++){
		vg->v[i] = v[i] + ref->v;
		vg->z[i] = z[i];
	}
	
	/* Update vg structure */
	for(i=n; i<(vg->n); i++)
		vg->v[i] = vg->z[i] = BADDOUBLE;
	vg->n = n;

on_error:
	
	/* Release memory */
	if(v) free(v);
	if(indx) free(indx);
	if(z) free(z);
	
	if(ierr)
		free_data((char *) vg, VGEO_VGEO_CODE, 1);

	return(ierr);

}
	
/* ------------------------------------------------------------------
	 -Main:
	 ------------------------------------------------------------------ */

#if PROTOTYPE_ALLOWED
void do_it(FILE *fp_cnt)
#else
void do_it()
#endif
{

	VGEO_TYPE         vgeo;
	VGEO_PRFT_TYPE    *prft;
	VGEO_CTD_TYPE     *ctd;
	VGEO_CTD_PRF_TYPE *prev, *curr;
	VGEO_VGEO_TYPE    *vg;
	VGEO_REF_TYPE     *ref;
	VGEO_FILES_TYPE   *fps;
	VGEO_VGRID_TYPE   *vgr;

	TIME_RANGE_TYPE trng;
	
	double v, e;
	int err, ierr=0;
	int i, n;

	char msg[200], str[200];

	/* Initializations */
	vgeo_init((VGEO_TYPE *) NULL);
	prft = &(vgeo.prft);
	ctd  = &(vgeo.ctd);
	prev = &(vgeo.ctd.prev);
	curr = &(vgeo.ctd.this);
	vg   = &(vgeo.vgeo);
	ref  = &(vgeo.ref);
	fps  = &(vgeo.fps);
	vgr  = &(vgeo.grid); 

	/* Get parameters */
	check_error(get_parameters(fp_cnt, vgeo_param, NULL),
							"Getting vgeo parameters");
	print_parameters(stdout, vgeo_param, NULL, "\n");
	
	/* Execute options */
	check_error(execute_options(fp_cnt, vgeo_options, ECHO) <= 0,
							"Reading vgeo options");

	/* Check options, initialize files and db */
	vgeo_init(&vgeo);

	/* Loop through time ranges */
	while(!ierr && 
				(get_time_range(fp_cnt, &(trng.start), &(trng.end)) == 1)) {
				
		/* Message */
		time_range_to_str(str, &(trng.start), &(trng.end));
		sprintf(msg, "\n%%%% TIME RANGE: %s\n", str);
		report_msg(msg);

		/* Sort CTD profiles for this time range */
		ierr = init_ctd_profiles(&trng, prft, year_base);
		ierr = db_not_eof(ierr);
		if(ierr && (ierr != EOF)) DBERROR(&ierr, "init_ctd_profiles()");
		if(!ierr || (ierr == EOF)){
			sprintf(msg, "%%%% Found %-d CTD profiles\n", prft->n);
			report_msg(msg);
		}
		/* Allow to skip to another time range */
		if(ierr == EOF){
			ierr = 0;
			continue;
		}

		/* Position ctd database inside time range and initialize first CTD
			 profile */
		prev->n = -1;
		while(!ierr && (prev->n <= 0)){
			/* Move to next CTD profile */
			ierr = next_ctd_profile(prft);
			err = db_not_eof(ierr);
			DBERROR(&err, "next_ctd_profile(): first");
			if(!ierr){
				/* Get current CTD profile */
				ierr = get_ctd_data(prev);
				if(ierr && (ierr != EOF)) DBERROR(&ierr, "get_ctd_data(): first");
				if(!ierr || (ierr == EOF)){
					yd_to_ymdhms_str(str, prev->yd, year_base);
					sprintf(msg, "%%%% CTD PROFLE TIME: %s", str);
					if(ierr) strcat(msg, " !!! NOT ENOUGH DATA !!!");
					strcat(msg, "\n");
					report_msg(msg);
					/* Allow to skip to next profile */
					if(ierr){
						ierr = 0;
						prev->n = -1;
						continue;
					}
				}
			}
		}
		/* Allow to skip to another time range */
		if(db_eof(ierr)){
			ierr = 0;
			continue;
		}
		
		/* Loop through the profiles */
		while(!ierr){
			
			/* Get next CTD data */
			curr->n = -1;
			while(!ierr && (curr->n <= 0)){
				/* Move to next CTD profile */
				ierr = next_ctd_profile(prft);
				err = db_not_eof(ierr);
				DBERROR(&err, "next_ctd_profile(): next");
				if(!ierr){
					/* Get current CTD profile */
					ierr = get_ctd_data(curr);
					if(ierr && (ierr != EOF)) DBERROR(&ierr, "get_ctd_data(): next");
					if(!ierr || (ierr == EOF)){
						yd_to_ymdhms_str(str, curr->yd, year_base);
						sprintf(msg, "%%%% CTD PROFLE TIME: %s", str);
						if(ierr) strcat(msg, " !!! NOT ENOUGH DATA !!!");
						strcat(msg, "\n");
						report_msg(msg);
						/* Allow to skip to next profile */
						if(ierr){
							ierr = 0;
							curr->n = -1;
							continue;
						}
					}
				}
			}
			
			if(!ierr){
				/* Calculate vertical derivative of geostrophic current */
				ierr = calculate_vgeo(vg, ctd);
				if(ierr && (ierr != EOF)) DBERROR(&ierr, "calculate_vgeo()");
				/* Allow to skip to next profile */
				if(ierr == EOF){
					yd_to_ymdhms_str(str, prev->yd, year_base);
					sprintf(msg, "%%%% VGEO: %s to ", str);
					yd_to_ymdhms_str(str, curr->yd, year_base);
					sprintf(msg, "%s: !!! NOT ENOUGH DATA !!!\n", str);
					report_msg(msg);
					ierr = copy_all_data((char *) prev, (char *) curr, VGEO_CTD_PRF_CODE);
					DBERROR(&ierr, "copy_all_data()");
					continue;
				}
			}

			if(!ierr){
				/* Get the reference level */
				ierr = get_reference(&vgeo);
				/* Allow to skip to next profile */
				if(db_eof(ierr)){
					yd_to_ymdhms_str(str, prev->yd, year_base);
					sprintf(msg, "%%%% VGEO: %s to ", str);
					report_msg(msg);
					yd_to_ymdhms_str(str, curr->yd, year_base);
					sprintf(msg, "%s: !!! NO REFERENCE !!!\n", str);
					report_msg(msg);
					ierr = copy_all_data((char *) prev, (char *) curr, 
															 VGEO_CTD_PRF_CODE);
					DBERROR(&ierr, "copy_all_data()");
					continue;
				}
				if(ierr) DBERROR(&ierr, "get_reference()\n");
			}

			if(!ierr){
				/* Integrate geostrophic current */
				ierr = integrate_vgeo(vg, ref);
				if(ierr && (ierr != EOF)) DBERROR(&ierr, "integrate_vgeo()");
				/* Allow to skip to next profile */
				if(ierr == EOF){
					ierr = copy_all_data((char *) prev, (char *) curr, VGEO_CTD_PRF_CODE);
					DBERROR(&ierr, "copy_all_data()");
					yd_to_ymdhms_str(str, prev->yd, year_base);
					sprintf(msg, "%%%% VGEO: %s to ", str);
					yd_to_ymdhms_str(str, curr->yd, year_base);
					sprintf(msg, "%s: !!! UNDEFINED REFERENCE !!!\n", str);
					report_msg(msg);
					continue;
				}
			}


			/* Print results */
			if(!ierr){
				if(CON_OUT){
					for(i=0; i<(vg->n); i++){
						v = BADDOUBLE;
						if(vg->v[i] < ADJ_BADFLOAT) v = vg->v[i]/opts.vgeo_scale;
						fprintf(fps->fp_con, "%15.8f %15.8f %15.8f %15.8f %15.8f\n",
										vg->yd, vg->pos.lon, vg->pos.lat, vg->z[i], v);
					}
				}
				if(STA_OUT){
					v = e = BADDOUBLE;
					if(ref->v < ADJ_BADFLOAT) v = ref->v/opts.vgeo_scale;
					if(ref->e < ADJ_BADFLOAT) e = ref->e/opts.vgeo_scale;
					fprintf(fps->fp_sta, 
									"%15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %6d\n",
									vg->yd, vg->pos.lon, vg->pos.lat, 
									ref->z, v, e, ref->n);
				}
				if(GRD_OUT){
					n = s_regrid_d(vg->z, vg->v, vgr->z, vgr->v, vg->n, vgr->n); 
					ierr = (n == EOF) ? INSUFFICIENT_MEMORY : 0;
					DBERROR(&ierr, "s_regrid_d(): regrid vg");
					if(!ierr){
						for(i=0; i<(vgr->n); i++){
							fprintf(fps->fp_grd, "%15.8f %15.8f %15.8f %15.8f",
											vg->yd, vg->pos.lon, vg->pos.lat, vgr->z[i]);
							if(vgr->v[i] < ADJ_BADFLOAT)
								fprintf(fps->fp_grd, " %15.8f\n", vgr->v[i]);
							else
								fprintf(fps->fp_grd, "          1.0E38\n");
						}
					}
				}
			}
			
			/* Move to next ensemble */
			if(!ierr){
				ierr = copy_all_data((char *) prev, (char *) curr, VGEO_CTD_PRF_CODE);
				DBERROR(&ierr, "copy_all_data()");
			}

		} /* Loop through CTD profiles */

		/* Check for errors */
		err = db_not_eof(ierr);
		ierr = DBERROR(&err, "Loop through CTD profiles");
		
				
	} /* Loop through CTD time ranges */
	
	/* Check for errors */
	err = db_not_eof(ierr);
	ierr = DBERROR(&err, "Loop through CTD time ranges");

	/* Exit status */
	if(!ierr)
		report_msg("%% VGEO terminated.\n");
	else
		DBERROR(&ierr, "Exiting VGEO");
	
	ierr = vgeo_exit(&vgeo);
	DBERROR(&ierr, "vgeo_exit()");

	return;

	
}

