/**--------------------------------------------------------------------------- ** ** 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 . If the table has an entry for the given , the corresponding options are executed. RETURNS: The value of . ------------------------------------------------------------------ */ #if PROTOTYPE_ALLOWED char *vgeo_read_options(/* File descriptor to the control file */ FILE *fp, /* Pointer to the table of correspondance between 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 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; iu = 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; iu) 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; iu[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; irho[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; iv[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; iz[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 , 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 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; itimes)+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 is defined, moves to the profile whose time is the first element of this list. If undefined, moves 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, , 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; ju)[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; iz[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; iunit_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; junit_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; ie)){ 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 and 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; iz[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; iv[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; }