/*************************************************************** ADCPSECT.C (derived from PROFSTAT.C) This program does lots of things with ADCP data in a CODAS database. It is still evolving, and I have not yet written a manual for it. ADCPSECT must be compiled with the LARGE MEMORY MODEL. Eric Firing Thu 11-03-1988 fixed bug in navigation option, so missing reference layer data is set to BADFLOAT. Sat 12-03-1988 Changed reference layer calculation to exclude data that are not within the good range as specified by the access variable structure (which counts bins starting from 1). Changed bin_range[] so that it counts bins from 1; subtract 1 when using as an index. Added time as an option for contouring x-coordinate. Made year_base a required parameter in the adcpsect parameter list, and removed it from the sequential_correlation and navigation option parameter lists. 89/01/18 Changed select_good() to mark velocities up to first_good_bin as BADFLOAT. Tue 02-07-1989 Broke the transport calculation out into a separate module. Sat 02-11-1989 Moved file name generation into a function, new_extension. Broke the sequential correlation function into a new module, seqcor.c; same for navigation and reference options, now in module shipref.c. Sun 02-12-1989 Moved contour file output to new module: confile.c. It still needs to be refined and reorganized. Sat 02-25-1989 More reorganization, including ctd.c; Added the option "nominal_svel". Added "integrate_top" to the regridding options, to cause the top bin value to be extrapolated to the surface, and all integration to start from the surface. Sat 03-11-1989 Added MATfile output to flux.c and to adcpsect.c, so that the ensemble mean times and positions will go in a file with extension ".mxy", the velocity will be in ".muv", and the transports in ".mtr". ".mxy": contains array "xyt" with a column for each ensemble, and three rows: longitude, latitude, decimal day; and vectors "z" and "zc" with the depth grid and the grid interval centers, respectively. ".muv": contains array "uv" with two columns for each ensemble, first u and then v. The row corresponds to the depth grid. ".mtr": contains array "transport" with two columns for each ensemble, first the integral of u dy, then -v dx. The vertical integral (cumulative summation) has not been done; do it in MATLAB with cumsum(transport). Note that the sign of the horizontal integral depends on the direction of ship motion. For example, with an eastward current the integral will be positive if the ship is going north, negative if southbound. Fri 03-17-1989 Fixed errors that had been made in adding the integrate_top option, which affected the gridding and editing regardless of which regrid option was used. Changed the handling of "output:" so that any extension is ignored, and the statistics output go to a file with extension ".sta". Wed 05-03-1989 Changed "write_navigation()" in SHIPREF.C so that all types of reference calculations can be used for .NAV output; previously "final_ship_ref" was omitted. In making the change I changed the algorithm for this function, making it much more compact. Tue 01-23-90 Changed time structure variable name from 'end' to (JR) 'end_time' because of problems on the Sun. Thu 08-16-90 Added PROFILE_FLAGS editing to select_good_velocity; (JR) Added flag_mask option to control file. Put call to select_good_velocity before call to get_ship.velocity Tue 09-24-91 For the rotate option, set ship heading = mean heading (JR) if available (nonzero for databases loaded after 9/13/91 correction to loadping), else use last heading Tue 07-28-92 Added SIGMA_T option to CTD module of adcpsect. (FB) Wed 09-22-93 Changed matlab output file extensions to _uv.mat, _xy.mat (JR) and _tr.mat because Matlab 4 assumes ASCII-formatted file for filenames with extension other than .mat Tue 08-09-94 Added "ascii" output option to derive ASCII subset with (PC) columns time, lon, lat, U & V at each depth layer Wed 12-07-94 Modified "ascii" output option to also include (PC) mean transducer temperature and ship velocity Fri 98/07/17 1) default for .sta ouput is none; "mean" must be specified explicitly. 2) New "time_ranges:" value "single" causes each individual profile to be written out. Tue 99/08/24 Added year_base to matlab output file. 2000/02/22 (EF) Added detide: option. This may be changed later for improved flexibility; at present it requires a specific ascii file format. Needed: Add "dir_position" code. Add Richardson number calculation, in module richno.c Add checking of data processing mask! ---------------------------------------------------------------- {control file structure:} {Comments and explanations are in curly brackets, and examples of control variables are in square brackets. Do not include the square brackets (or curly ones) in the non-commented (active) part of an actual control file. C-style comments can be used anywhere, so long as they are not nested. Anything that is not in either kind of brackets is a key word.} dbname: [WEP] output: [wepsect] {no extension; all extensions will be added automatically} step_size: [1] {number of profiles to advance} ndepth: [60] {number of depth bins to read and process} time_ranges: [combined] { or separate or single ) year_base= [1988] option_list: verbose { prints lots of debugging and auxiliary output } underway? [ yes ] limit= [ 3 ] { limit in m/s } pg_min= [ 80 ] AGC_margin= [ 7 ] skew_limit: agc_limit= [100] max_velocity_difference= [ 0.6 ] {m/s} reference: [reference_bins] [ 5 ] to [ 20 ] { bin range; first bin is 1 } { or one of the following: none final_ship_ref nav_position nav_velocity bottom_track dir_position } {dir_position not available yet} navigation: [reference_bins] [ 5 ] to [ 20 ] { bin range; first bin is 1 } { and/or any of the following: final_ship_ref nav_position nav_velocity bottom_track dir_position } {dir_position not available yet} end {required to end navigation option list} detide: filename: [ tide.asc ] { ascii file with columns: } { decimal day } { decimal longitude } { decimal latitude } { predicted tide U, mm/s} { predicted tide V, mm/s} { predicted tide h, mm } rotate: [ 45 ] { angle in degrees, old x-axis to new x-axis } { directions are counterclockwise } { or substitute "ship_coordinates"} ctd: dbname: [ \wepocs\ctd\c2us ] nominal_svel= [1470.0] {default, for RDI; use 1500 for Ametek} regrid: [integrate ] { or integrate_top, interpolate, centered_average, average, per_grid_interval } [ depth ] { or svan, pot_svan, or sigma_theta (if ctd option selected) } {must have one and only one of the following: } grid number= [ 50 ] origin= [ 16 ] increment= [ 8 ] grid_list number= [ 5 ] boundaries: [20] [40] [60] [80] [100] contour: [ latitude ] { or longitude or time} [ middle ] { or mean } minimum_npts= [ 10 ] units= [0.01] vector: minimum_npts= [ 10 ] ascii: { time, position, u v repeated } minimum_percent= [ 0.50 ] { %ensembles in range with good data } { less than this, program prints bad = 99999 } end {required to end ascii option list} transport: minimum_npts= [ 10 ] units= [1.e6] {for Sverdrups} direction= [all] {or left or right} sequential_cor: bin_range: [3] to [40] {starting from 1, inclusive} statistics: { followed by list of options, terminated by end } [mean] { and/or any or all of the following: variance, extrema, norm_extrema, indices, all} units= [0.01] end {required to end statistics option list} flag_mask: [ bit-mask ] ... {see adcp/include/adcpflag.h} end {required to end flag mask list} end {required to end adcpsect option list} {list of YMDHMS time pairs:} {(yy/mm/dd hh:mm:ss) to (yy/mm/dd hh:mm:ss)} {etc.} ****************************************************************/ #include "geninc.h" #include #include "ioserv.h" #include "vector.h" #include "use_db.h" #include "vstat.h" #include "ocean.h" #include "pos_to_m.h" #include "adcp.h" #include "adcpflag.h" #include "adcpsect.h" #include "data_id.h" #ifndef M_PI #define M_PI 3.14159265358979323846 #endif /*--> ADCP general parameters */ YMDHMS_TIME_TYPE start, end_time; /* for time range */ int verbose = FALSE; /* flag for debugging output to stdout */ int first_profile; /* flag: true if first profile in ensemble */ char time_stamp[30]; /* global variable set to start time of * calculation. */ /* It will be set in do_it, using ctime(time). */ static FILE_NAME_TYPE DBNAME; static int STEPS; /* used by write_ensemble() */ static int ndepth; /* max number of depth bins to process */ static int combined = 0; /* flag: 1 if all time ranges in one ensemble */ static int separate = 0; /* flag: 1 if one time range per ensemble */ static int single_ens = 0; /* flag: 1 if one profile per ensemble */ /* (Note: ident "single" conflicts with a Solaris header file declaration! */ static int rotate = FALSE; /* flag */ static int on_station = FALSE; /* flag: select on-station profiles */ extern int use_ctd; /* flag; ctd data in use */ extern int detide; /* flag; detiding is active */ static int statistics = FALSE; /* former default was U_MEAN; multi-bit flag */ static int vector = FALSE; /* flag */ extern int transport; /* flag in flux.c */ extern int seq_cor; /* flag in seqcor.c */ extern int contour; /* flag */ static int ascii; /* flag */ static ADCP_PROFILE_TYPE adcp; static CTD_PROFILE_TYPE ctd; /*--> ADCP statistics */ static UNISTAT_TYPE u_stat, v_stat; static double stat_units = 1.0; UNISTAT_TYPE adcp_pos_stat; /* longitude, latitude, decimal day */ UNISTAT_TYPE adcp_ship_stat; /* ship's velocity u and v */ UNISTAT_TYPE tr_temp_stat; /* transducer temperatute */ double mid_adcp_pos[3]; /* avg of min and max */ static double *temp_d_array; /* for conversion to NaN flagging */ /*--> ASCII standard subset file */ static float asc_min_pct = 0.50;/* minimum number for ascii output */ NAME_LIST_ENTRY_TYPE lat_lon_names[] = { {"longitude", LONGITUDE}, {"latitude", LATITUDE}, {"time", DDTIME}, {NULL, 0} }; static OPTION_TYPE unistat_options[] = { {"end", 0, NULL, NULL}, {"mean", U_MEAN, logical_or, (char *) &statistics}, {"variance", U_VARIANCE, logical_or, (char *) &statistics}, {"extrema", U_EXTREMA, logical_or, (char *) &statistics}, {"norm_extrema", U_NORM_EXTREMA, logical_or, (char *) &statistics}, {"indices", U_INDICES, logical_or, (char *) &statistics}, {"all", U_ALL, logical_or, (char *) &statistics}, {"units=", TYPE_DOUBLE, op_get_param, (char *) &stat_units}, {NULL, 0, NULL, NULL} }; static OPTION_TYPE ascii_options[] = { {"end", 0, NULL, NULL}, {"minimum_percent=", TYPE_FLOAT, op_get_param, (char *) &asc_min_pct}, {NULL, 0, NULL, NULL} }; /*--> rotation parameters */ static double angle_old_to_new; /* angle between old and new coords */ #define ROTATE_ABSOLUTE 1 #define ROTATE_SHIP 2 /*--> parameters for profile selection and editing */ static double speed_limit = 0.0;/* criterion for on or off station */ static int pg_min = 0; /* threshold for data acceptance */ static int amp_margin = 0; /* threshold for data acceptance */ static int agc_limit = 0; static UBYTE flag_mask = ALL_BITS; /* PROFILE_FLAGS bits to consider jr+ */ static double bad_dif = 0.0; /*--> regridding variables and constants: */ extern float sigma_sva_grid[NDMAX_NG]; extern float adj_depth[NDMAX_NG]; extern float adj_u[NDMAX_NG]; extern float adj_v[NDMAX_NG]; extern int adj_ndepth; /* new grid; = rg_number if regridding is * done */ extern int adj_nu; /* adj_ndepth or one less; npts in array of * centers */ extern VERTICAL_GRID_TYPE vgrid; extern int regrid_var; extern int regrid_method; extern NAME_LIST_ENTRY_TYPE regrid_methods[]; extern NAME_LIST_ENTRY_TYPE regrid_variables[]; /* * parameter in ctd.c used for depth correction. It is global so that it can * be set easily by an option in the main option list. It would more * logically be local to ctd.c, and part of a local option list there, but * that would require a change to the control file format. */ extern double nominal_svel; static int year_base; extern PROFILE_REF_TYPE ref; extern PROFILE_REF_TYPE nav; /*--> output parameters and variables */ FILE_NAME_TYPE out_filename; /* global, so other modules can use it to * construct file names. In the other * modules, it must be declared as extern * char out_filename[]; declaring it as a * pointer does not work. */ static FILE_NAME_TYPE vec_filename, stat_filename, mat_filename, asc_filename; static FILE *fp_out, *fp_zxyt_mat, *fp_uv_mat, *fp_vec, *fp_asc; static int vec_min_npts; /* minimum number for vector output */ static MAT_ARRAY_FP_TYPE *mat_zxyt_ptr, *mat_uv_ptr; /*************************************************************** Parameter input structures for do_it in ADCPSECT: */ NAME_LIST_ENTRY_TYPE time_ranges[] = { {"separate", FALSE}, {"combined", TRUE}, {"single", -1}, {NULL, 0} }; NAME_LIST_TYPE adcpsect_name_lists[] = { time_ranges }; static PARAMETER_LIST_ENTRY_TYPE adcpsect_params[] = { {" dbname:", " %79s", DBNAME, TYPE_STRING}, {" output:", " %79s", out_filename, TYPE_STRING}, {" step_size:", " %d", &STEPS, TYPE_INT}, {" ndepth:", " %d", &ndepth, TYPE_INT}, {" time_ranges:", " %19s", &combined, TYPE_TRANS}, {" year_base=", " %d", &year_base, TYPE_INT}, {" option_list:", NULL, NULL, 0}, {NULL, NULL, NULL, 0} }; #if PROTOTYPE_ALLOWED int speed_in_range(ADCP_PROFILE_TYPE * p, double ship_u, double ship_v) #else int speed_in_range(p, ship_u, ship_v) ADCP_PROFILE_TYPE *p; double ship_u, ship_v; #endif /* returns TRUE if ship speed is in range */ { double speed; if (speed_limit == 0.0) return (TRUE); /* no speed checking */ if (good_float(ship_u)) /* a ship velocity is available; use it */ speed = sqrt(sqr(ship_u) + sqr(ship_v)); else if (good_float(p->u[2]))/* otherwise use relative velocity in 3rd bin */ speed = sqrt(sqr(p->u[2]) + sqr(p->v[2])); else return (FALSE); /* profile is not much good anyway */ if (on_station) /* select only on-station profiles */ { if (speed < speed_limit) return (TRUE); else return (FALSE); } else /* select underway profiles */ { if (speed >= speed_limit) return (TRUE); else return (FALSE); } } /* speed_in_range */ #if PROTOTYPE_ALLOWED int index_bad_tracking(ADCP_PROFILE_TYPE * p, int agc_limit, double bad_dif) #else int index_bad_tracking(p, agc_limit, bad_dif) ADCP_PROFILE_TYPE *p; int agc_limit; double bad_dif; #endif { int i; int itop, n0; struct { double u, v; } pmin, pmax, pdif; double magdif = 0.0; /* magnitude of difference */ pmin.u = pmin.v = BADFLOAT; pmax.u = pmax.v = -BADFLOAT; n0 = min_val(p->namp, p->nu); for (itop = 1; itop < n0 - 1 && p->amp[itop + 1] > agc_limit; itop++); /* ef */ for (i = itop; i < n0 && magdif <= bad_dif; i++) { if (p->u[i] < ADJ_BADFLOAT) { pmin.u = min_val(pmin.u, p->u[i]); pmin.v = min_val(pmin.v, p->v[i]); pmax.u = max_val(pmax.u, p->u[i]); pmax.v = max_val(pmax.v, p->v[i]); pdif.u = pmax.u - pmin.u; pdif.v = pmax.v - pmin.v; magdif = sqrt((pdif.u * pdif.u) + (pdif.v * pdif.v)); } } if (i == n0) return (0); /* threshold was not reached */ else return (i); /* first index where threshold was met or * exceeded */ } /* index_bad_tracking() */ #if PROTOTYPE_ALLOWED void select_good_velocity(ADCP_PROFILE_TYPE * p) #else void select_good_velocity(p) ADCP_PROFILE_TYPE *p; #endif /* this can be expanded to use various criteria */ { int i; unsigned int npf; /* jr+ */ p->namp = -1; /* flag to show it has not been read *//* ef */ if (p->access.last_good_bin < MAXSHORT) /* check added 97/04/25 */ p->nu = min_val(p->nu, p->access.last_good_bin + 1); /* ef */ /* This check is needed only because for a time, updscan and dbupdate could set last_good_bin to MAXSHORT; they have now been changed to set it to its default, conf1.num_bins, instead. */ for (i = 1; i < (p->access.first_good_bin); i++) /* ef */ { p->u[i] = BADFLOAT; p->v[i] = BADFLOAT; } if (pg_min) { p->npg = 1 + check_dbget(PERCENT_GOOD, p->percent_good + 1, ndepth, /* ef */ "DBGET, percent good"); p->nu = min_val(p->nu, p->npg); /* * Reduce nu until all bad points at the end of the profile are * eliminated. The minimum number of points is 1, u[0], which is not * read in from the data but may be used to extrapolate to the surface * (by setting it equal to the top bin). */ while ((p->percent_good[p->nu - 1] < pg_min) && p->nu > 1) p->nu--; /* Now flag any bad points in the middle of the profile. */ for (i = 1; i < p->nu; i++) /* ef */ { if (p->percent_good[i] < pg_min) { p->u[i] = BADFLOAT; p->v[i] = BADFLOAT; } } } /* end pg_min screening */ if (amp_margin) { int noise_floor = 0; int count = 0; int min_amp; p->namp = 1 + check_dbget(AMP_SOUND_SCAT, (char *) (p->amp + 1), ndepth, "DBGET, amplitude"); /* Define the noise floor as average AGC over last 10 bins. */ for (i = p->namp - 10; i < p->namp; i++) /* ef */ { if (p->amp[i] != BADUBYTE) { noise_floor += p->amp[i]; count++; } } if (count) noise_floor /= count; min_amp = noise_floor + amp_margin; p->nu = min_val(p->nu, p->namp); while ((p->amp[p->nu - 1] < min_amp) && p->nu > 1) p->nu--; for (i = 1; i < p->nu; i++) /* ef */ { if (p->amp[i] < min_amp) { p->u[i] = BADFLOAT; p->v[i] = BADFLOAT; } } } /* end amp_margin screening */ if (agc_limit) { int ibad; if (p->namp < 0) /* ef */ p->namp = 1 + check_dbget(AMP_SOUND_SCAT, (char *) (p->amp + 1), ndepth, "DBGET, amplitude"); ibad = index_bad_tracking(p, agc_limit, bad_dif); if (ibad) { p->nu = ibad; printf("\nskew error starting at %d ", ibad); } } /* end skew error screening */ if (flag_mask && p->nu > 0 &&/* jr+ */ (0 < (npf = check_dbget(PROFILE_FLAGS, (char *) (p->profile_flags + 1), (p->nu), "DBGET, profile flags")))) /* use p->nu to * immediately limit npf */ { flag_data(p->u + 1, p->profile_flags + 1, npf, flag_mask); flag_data(p->v + 1, p->profile_flags + 1, npf, flag_mask); } /* end profile_flags screening */ } /* select_good_velocity */ #if PROTOTYPE_ALLOWED void initialize_ensemble(void) #else void initialize_ensemble() #endif { zero_unistat(&u_stat); zero_unistat(&v_stat); zero_unistat(&adcp_pos_stat); zero_unistat(&adcp_ship_stat); zero_unistat(&tr_temp_stat); adcp.depth[0] = 0.0; first_profile = TRUE; } /* initialize_ensemble */ /* * 91/09/24 JR- : use mean_heading if available (loaded after 91/09/13 * correction to loadping) else use last_heading #define SHIP_HEADING * adcp.ancillary_2.last_heading */ #if PROTOTYPE_ALLOWED void increment_ensemble(void) #else void increment_ensemble() #endif { int i; int nu, nv; float *u_ptr, *v_ptr; static FIX_TYPE last_fix; struct { double u, v; } ship; double ship_heading; check_dbget(ANCILLARY_1, (char *) &adcp.ancillary_1, sizeof(ANCILLARY_1_TYPE), "DBGET, ancillary_1"); check_dbget(ANCILLARY_2, (char *) &adcp.ancillary_2, sizeof(ANCILLARY_2_TYPE), "DBGET, ancillary_2"); check_dbget(ACCESS_VARIABLES, (char *) &adcp.access, sizeof(ACCESS_VARIABLES_TYPE), "DBGET, access"); /* * We have to get u and v arrays in order to get the ship velocity and * check for speed range (in case navigation is not available, or the bin * reference option has been chosen). */ nu = check_dbget_f(U, adcp.u + 1, ndepth, "DBGET_F, u"); nv = check_dbget_f(V, adcp.v + 1, ndepth, "DBGET_F, v"); adcp.nu = 1 + min_val(nu, nv); adcp.u[0] = BADFLOAT; adcp.v[0] = BADFLOAT; if (detide) subtract_tide(&adcp); get_ship_velocity(&adcp, &ref, &last_fix, &ship.u, &ship.v); if (verbose) printf(" ship U= %f V= %f; ", ship.u, ship.v); if (speed_in_range(&adcp, ship.u, ship.v)) /* otherwise,skip everything */ { adcp.nd = 1 + check_dbget_f(DEPTH, adcp.depth + 1, /* ef */ ndepth, "DBGET_F, depth"); if (use_ctd) /* We should also check the data processing * mask. */ correct_depths(&ctd, adcp.depth, &adcp.nd); if (adcp.nd < adcp.nu) adcp.nu = adcp.nd; /* 94/02/10 jr+ */ select_good_velocity(&adcp); if (nav.method != REF_NONE) write_navigation(&adcp, &nav); get_latlon(&adcp.position); update_unistat_d(&adcp_pos_stat, (double *) &adcp.position, 3); update_unistat_d(&adcp_ship_stat, (double *) &ship, 2); update_unistat(&tr_temp_stat, &adcp.ancillary_1.tr_temp, 1); if (ref.method != REF_NONE) /* we should check data proc mask also */ { if (good_float(ship.u)) { for (i = 0, u_ptr = adcp.u; i < adcp.nu; i++, u_ptr++) if (*u_ptr < ADJ_BADFLOAT) *u_ptr += ship.u; for (i = 0, v_ptr = adcp.v; i < adcp.nu; i++, v_ptr++) if (*v_ptr < ADJ_BADFLOAT) *v_ptr += ship.v; } else { for (i = 0, u_ptr = adcp.u; i < adcp.nu; i++, u_ptr++) *u_ptr = BADFLOAT; for (i = 0, v_ptr = adcp.v; i < adcp.nu; i++, v_ptr++) *v_ptr = BADFLOAT; } } if (rotate) { if (rotate == ROTATE_SHIP) /* v is forward, u is starboard */ { if (adcp.ancillary_1.mn_heading == 0.0) ship_heading = (double) adcp.ancillary_2.last_heading; else ship_heading = (double) adcp.ancillary_1.mn_heading; angle_old_to_new = -(M_PI / 180) * ship_heading; } if (verbose) printf(" rotation angle: %7.2f radians ", angle_old_to_new); rotate_vector(angle_old_to_new, adcp.u, adcp.v, adcp.nu); } regrid_velocity(&adcp, &ctd); if (seq_cor) calc_seqcor(adj_u, adj_v, adj_nu, adcp.ddtime); update_unistat(&u_stat, adj_u, adj_nu); update_unistat(&v_stat, adj_v, adj_nu); first_profile = FALSE; } } /* increment_ensemble */ #if PROTOTYPE_ALLOWED void calculate_ensemble(void) #else void calculate_ensemble() #endif { int i; calculate_unistat(&u_stat); calculate_unistat(&v_stat); calculate_unistat(&adcp_pos_stat); for (i = 0; i < 3; i++) { mid_adcp_pos[i] = adcp_pos_stat.npts[i] > 0 ? 0.5 * (adcp_pos_stat.min[i] + adcp_pos_stat.max[i]) : BADFLOAT; } if (transport) calculate_transport(&u_stat, &v_stat, &adcp_pos_stat); if (ascii) { calculate_unistat(&adcp_ship_stat); calculate_unistat(&tr_temp_stat); if (tr_temp_stat.sum[0] <= ADJ_BADFLOAT) if (tr_temp_stat.sum[0] < -10. || tr_temp_stat.sum[0] > 37.) { printf("\n\n EXTREME TEMPERATURE = %8.2f time=%9.5f\n", tr_temp_stat.sum[0], adcp_pos_stat.sum[DDTIME]); /* tr_temp_stat.sum[0] = ADJ_BADFLOAT; */ } if (verbose) { printf(" \nship_stat=%f sd=%f\n", adcp_ship_stat.sum[0], adcp_ship_stat.sumsq[0]); printf("\n tr=%8.2f sd=%8.2f\n", tr_temp_stat.sum[0], tr_temp_stat.sumsq[0]); } } } /* calculate_ensemble */ void end_ensemble(void) { if (verbose) printf("\nready to calculate\n"); calculate_ensemble(); if (verbose) printf("\nready to print\n"); print_ensemble(); } #if PROTOTYPE_ALLOWED void print_header(FILE * fp) #else void print_header(fp) FILE *fp; #endif { fprintf(fp, "\nADCPSECT %s\n", time_stamp); print_parameters(fp, adcpsect_params, adcpsect_name_lists, "\n"); if (use_ctd) write_ctd_params(fp); write_vgrid_params(fp); if (ref.method == REF_BINS) fprintf(fp, "\nReference layer: bins %d through %d", ref.bin_range[0], ref.bin_range[1]); else fprintf(fp, "\nNo reference layer removed."); if (speed_limit != 0.0) fprintf(fp, "\nProfiles from %s, speed limit = %f m/s.", (on_station ? "on station" : "underway"), speed_limit); if (!pg_min) fprintf(fp, "\n minimum percent good = %d", pg_min); if (!amp_margin) fprintf(fp, "\n minimum AGC counts above noise floor = %d", amp_margin); if (rotate == ROTATE_ABSOLUTE) fprintf(fp, "\n x-axis rotated %f degrees CCW from east.", angle_old_to_new * 180.0 / M_PI); else if (rotate == ROTATE_SHIP) fprintf(fp, "\n Rotation: x-axis is starboard, y is forward."); if (detide) write_detide_params(fp); } /* print_header() */ #define print_fixed(x,w,p) x>=ADJ_BADFLOAT ? fprintf(fp_vec, " %*s", w, "1E38") : fprintf(fp_vec, " %*.*f", w, p, x) #if PROTOTYPE_ALLOWED void print_ensemble(void) #else void print_ensemble() #endif { int i; int min_matuv = 1; /* This could be made a parameter. */ if (verbose) printf("\nentering print_ensemble\n"); if (contour) write_confile(&u_stat, &v_stat, vgrid.center_array, vgrid.increase_down); if (vector) /* * This is a minimal section to put out a file that can be input to a * program to plot velocity (or transport) arrows. The main ASCII file * will have to be used as documentation of the parameters. If * adj_ndepth is not small (like a half-dozen), the lines will be very * long! The format for each line is: * * longitude latitude u1 v1 u2 v2 ....... * * MKS units are used. */ { if (adcp_pos_stat.sum[LONGITUDE] <= ADJ_BADFLOAT) fprintf(fp_vec, "%9.4f %9.4f", adcp_pos_stat.sum[LONGITUDE], adcp_pos_stat.sum[LATITUDE]); else fprintf(fp_vec, " 1E38 1E38"); for (i = 0; i < adj_nu; i++) { if (u_stat.npts[i] >= vec_min_npts) { fprintf(fp_vec, " %6.2f", u_stat.sum[i]); fprintf(fp_vec, " %6.2f", v_stat.sum[i]); } else { fprintf(fp_vec, " 1E38 1E38"); } } fprintf(fp_vec, "\n"); } if (ascii) /* * MKS units are used. */ { /* time as decimal day */ if (adcp_pos_stat.sum[DDTIME] <= ADJ_BADFLOAT) fprintf(fp_asc, "%9.5f", adcp_pos_stat.sum[DDTIME]); else /* use middle of input time range as * placeholder for missing data */ fprintf(fp_asc, "%9.5f", ((year_day(&start, year_base) + year_day(&end_time, year_base)) / 2.)); /* write out position as longitude and latitude; missing = 1E38 */ if (adcp_pos_stat.sum[LONGITUDE] <= ADJ_BADFLOAT) { fprintf(fp_asc, " %9.4f %8.4f", (adcp_pos_stat.sum[LONGITUDE] <= 180 ? adcp_pos_stat.sum[LONGITUDE] : adcp_pos_stat.sum[LONGITUDE] - 360.0), adcp_pos_stat.sum[LATITUDE]); } else fprintf(fp_asc, " 1E38 1E38"); /* * write out transducer temperature (deg C) mean and variance missing = * 1E38 */ if (tr_temp_stat.sum[0] <= ADJ_BADFLOAT) fprintf(fp_asc, " %4.1f %5.2f", tr_temp_stat.sum[0], sqrt(tr_temp_stat.sumsq[0])); else fprintf(fp_asc, " 1E38 1E38"); /* write out ship mean speed and variance; missing = 1E38 (m/sec) */ if (adcp_ship_stat.sum[0] <= ADJ_BADFLOAT) fprintf(fp_asc, " %5.1f %5.2f %5.1f %5.2f", adcp_ship_stat.sum[0], sqrt(adcp_ship_stat.sumsq[0]), adcp_ship_stat.sum[1], sqrt(adcp_ship_stat.sumsq[1])); else fprintf(fp_asc, " 1E38 1E38 1E38 1E38"); /* write out each depth pair of u and v components */ for (i = 0; i < adj_nu; i++) { if ( ((float) u_stat.npts[i]) / u_stat.n_increments >= asc_min_pct && u_stat.sum[i] < ADJ_BADFLOAT) { fprintf(fp_asc, " %5.0f %5.0f", u_stat.sum[i] / 0.001, v_stat.sum[i] / 0.001); } else { fprintf(fp_asc, " 99999 99999"); } } fprintf(fp_asc, "\n"); } /* if (ascii) */ if (statistics) { /* always print the time and position information */ if (separate) print_time_range(fp_out, &start, &end_time); fprintf(fp_out, "\n ADCP mean position: %g lat, %g long", adcp_pos_stat.sum[LATITUDE], adcp_pos_stat.sum[LONGITUDE]); fprintf(fp_out, "\n ADCP middle position: %g lat, %g long, %g days", mid_adcp_pos[LATITUDE], mid_adcp_pos[LONGITUDE], mid_adcp_pos[DDTIME]); fprintf(fp_out, "\n ADCP latitude range: %g to %g", adcp_pos_stat.min[LATITUDE], adcp_pos_stat.max[LATITUDE]); fprintf(fp_out, "\n ADCP longitude range: %g to %g", adcp_pos_stat.min[LONGITUDE], adcp_pos_stat.max[LONGITUDE]); if (use_ctd) { fprintf(fp_out, "\n CTD position: %f %f time: ", ctd.position.lat, ctd.position.lon); write_ymdhms_time(fp_out, &(ctd.time)); } fprintf(fp_out, "\n No. of ensembles in time range: %d\n", u_stat.n_increments); write_unistat(fp_out, &u_stat, stat_units, vgrid.array, vgrid.label, vgrid.precision); write_unistat(fp_out, &v_stat, stat_units, vgrid.array, vgrid.label, vgrid.precision); fprintf(fp_out, "\n\n"); } /* if (statistics) */ if (transport) print_transport(&adcp_pos_stat, vgrid.array, vgrid.label, vgrid.precision); tonan_d(temp_d_array, adcp_pos_stat.sum, 3); add_to_mat_array(mat_zxyt_ptr, (char *) temp_d_array); tonan_mask_d(temp_d_array, u_stat.sum, u_stat.npts, min_matuv, adj_nu); add_to_mat_array(mat_uv_ptr, (char *) temp_d_array); tonan_mask_d(temp_d_array, v_stat.sum, v_stat.npts, min_matuv, adj_nu); add_to_mat_array(mat_uv_ptr, (char *) temp_d_array); } /* print_ensemble */ /*************************************************************** Option initialization functions: */ #if PROTOTYPE_ALLOWED int init_stats(FILE * fp_cnt) #else int init_stats(fp_cnt) FILE *fp_cnt; #endif { check_error((execute_options(fp_cnt, unistat_options, ECHO) <= 0), "reading statistics options"); return (0); } #if PROTOTYPE_ALLOWED int init_uw(FILE * fp_cnt) #else int init_uw(fp_cnt) FILE *fp_cnt; #endif { char sub_option[5]; check_error((fscanf(fp_cnt, " %4s limit= %lf", sub_option, &speed_limit) != 2), "reading speed limit"); if (strcmp(sub_option, "yes") == 0) on_station = FALSE; else if (strcmp(sub_option, "no") == 0) on_station = TRUE; else check_error(1, "underway? yes or no"); return (0); } #if PROTOTYPE_ALLOWED int init_rotate(FILE * fp_cnt) #else int init_rotate(fp_cnt) FILE *fp_cnt; #endif { char sub_option[20]; check_error((fscanf(fp_cnt, " %19s", sub_option) != 1), "reading rotation option"); if (!strcmp(sub_option, "ship_coordinates")) { rotate = ROTATE_SHIP; } else /* should be absolute rotation angle */ { check_error( (sscanf(sub_option, " %lf", &angle_old_to_new) != 1), "reading rotation angle"); angle_old_to_new *= M_PI/180.0; rotate = ROTATE_ABSOLUTE; } return(0); } #if PROTOTYPE_ALLOWED int init_vector(FILE *fp_cnt) #else int init_vector(fp_cnt) FILE *fp_cnt; #endif { check_error( (fscanf(fp_cnt, " minimum_npts= %d", &vec_min_npts) != 1), "reading vector parameter"); vector = TRUE; /* Make a vector output file name, and open the file. */ new_extension(vec_filename, out_filename, ".vec"); fp_vec = check_fopen(vec_filename, "w"); return (0); } #if PROTOTYPE_ALLOWED int init_ascii(FILE * fp_cnt) #else int init_ascii(fp_cnt) FILE *fp_cnt; #endif { int i; check_error((execute_options(fp_cnt, ascii_options, ECHO) <= 0), "reading ascii options"); ascii = TRUE; /* Make a ascii output file name, and open the file. */ new_extension(asc_filename, out_filename, ".sub"); fp_asc = check_fopen(asc_filename, "w"); /* * Scan each block file to find maximum number of levels and minimum start * level */ /* Make header */ for (i = 1; i < NDMAX_NG; i++) if (adj_depth[i] == 0.0) break; /* 2002/01/17 Pierre Jaccard * Prepend a '%' to the header */ if (ref.method == REF_FINAL_SHIP) fprintf(fp_asc, "%%sac_id=%s yr_base=%d start_lev=%3.0fm num_lev=%3d absolute\n", out_filename, year_base, adj_depth[0], i); else fprintf(fp_asc, "%%sac_id=%s yr_base=%d start_lev=%3.0fm num_lev=%3d relative\n", out_filename, year_base, adj_depth[0], i); return (0); } #if PROTOTYPE_ALLOWED int init_skew(FILE * fp_cnt) #else int init_skew(fp_cnt) FILE *fp_cnt; #endif { check_error((fscanf(fp_cnt, " agc_limit= %d max_velocity_difference= %lf", &agc_limit, &bad_dif) != 2), "reading skew parameters"); return (0); } #if PROTOTYPE_ALLOWED int init_mask(FILE * fp_cnt) /* jr+ */ #else int init_mask(fp_cnt) /* jr+ */ FILE *fp_cnt; #endif { char flag_name[40]; int bit_flag; flag_mask = 0; do { check_error(!fscanf(fp_cnt, " %s", flag_name), "reading flag_mask"); if (!strcmp(flag_name, "end")) break; check_error((BADINT == (bit_flag = get_code(adcp_profflags, flag_name))), "unrecognized flag_mask"); flag_mask |= (UBYTE) bit_flag; /* cumulate flag bit-by-bit */ } while (1); return (0); } static OPTION_TYPE adcpsect_options[] = { {"end", 0, NULL, (char *) NULL}, {"reference:", 1, init_ref, (char *) NULL}, {"navigation:", 2, init_nav, (char *) NULL}, {"pg_min=", TYPE_INT, op_get_param, (char *) &pg_min}, {"verbose", 4, set_code, (char *) &verbose}, {"underway?", 5, init_uw, (char *) NULL}, {"rotate:", 6, init_rotate, (char *) NULL}, {"ctd:", 7, init_ctd, (char *) NULL}, {"nominal_svel=", TYPE_DOUBLE, op_get_param, (char *) &nominal_svel}, {"regrid:", 8, init_regrid, (char *) NULL}, {"statistics:", 9, init_stats, (char *) NULL}, {"contour:", 10, init_contour, (char *) NULL}, {"vector:", 11, init_vector, (char *) NULL}, {"ascii:", 16, init_ascii, (char *) NULL}, {"transport:", 13, init_transport, (char *) NULL}, {"sequential_cor:", 14, init_seqcor, (char *) NULL}, {"AGC_margin=", TYPE_INT, op_get_param, (char *) &_margin}, {"skew_limit:", 12, init_skew, (char *) NULL}, {"flag_mask:", 15, init_mask, (char *) NULL}, {"detide:", 17, init_detide, (char *) NULL}, {NULL, 0, NULL, NULL} }; /*************************************************************** Main routine: read the control file, process the data accordingly. */ #if PROTOTYPE_ALLOWED void do_it(FILE * fp_cnt) #else void do_it(fp_cnt) FILE *fp_cnt; #endif { int IERR = 0, nbin, i; time_t clock_start, clock_end; int run_time; /* extern long timezone; *//* C library global variables */ /* extern int daylight; -- not needed & machine-dependent */ /* timezone = 10 * 60 * 60; *//* for Hawaii */ /* daylight = 0; */ time(&clock_start); strncpy(time_stamp, ctime(&clock_start), 24); /* stop at 24 to strip * \n */ check_error(get_parameters(fp_cnt, adcpsect_params, adcpsect_name_lists), "getting adcpsect parameters"); if (combined == -1) { single_ens = 1; combined = 0;} else if (combined == FALSE) {separate = 1; } check_dbopen(ADCP_DB, DBNAME, READ_ONLY, DIR_IN_MEMORY); printf("\n\n DATABASE: %s SUCCESSFULLY OPENED", DBNAME); new_extension(stat_filename, out_filename, ".sta"); fp_out = check_fopen(stat_filename, "w"); new_extension(mat_filename, out_filename, "_xy.mat"); fp_zxyt_mat = check_fopen(mat_filename, "wb"); new_extension(mat_filename, out_filename, "_uv.mat"); fp_uv_mat = check_fopen(mat_filename, "wb"); /* * ---> Read the list of options. */ if (execute_options(fp_cnt, adcpsect_options, ECHO) <= 0) { printf("\n ERROR in reading adcpsect options."); exit(-1); } check_error((regrid_var != RV_DEPTH && regrid_var != RV_NONE) && !use_ctd, "regridding by svan or sigma theta requires ctd option"); /* * ---> Calculate array sizes and allocate the arrays. */ ndepth = min_val(ndepth, NDMAX - 1); /* ef */ set_vgrid(ndepth); /* sets adj_ndepth and adj_nu */ set_ctd_grid(&ctd, regrid_var); allocate_unistat(&u_stat, adj_nu, "U", statistics); allocate_unistat(&v_stat, adj_nu, "V", statistics); allocate_unistat(&adcp_pos_stat, 3, "ADCP positions", U_MEAN | U_EXTREMA); allocate_unistat(&tr_temp_stat, 1, "TR-TEMP", U_MEAN | U_VARIANCE); allocate_unistat(&adcp_ship_stat, 2, "SHIP SPEED", U_MEAN | U_VARIANCE); if (transport) allocate_transport(adj_nu); nbin = max_val(3, adj_ndepth); temp_d_array = (double *) calloc(nbin, sizeof(double)); check_error((temp_d_array == NULL), "insufficient memory for temp_d_array"); /* * 91/11/18 jr++: copy all non-double to double before saving to * mat-file: */ for (i = 0; i < adj_ndepth; i++) temp_d_array[i] = (double) vgrid.array[i]; savemat(fp_zxyt_mat, 'd', "z", adj_ndepth, 1, 0, (char *) temp_d_array, (char *) NULL); for (i = 0; i < adj_nu; i++) temp_d_array[i] = (double) vgrid.center_array[i]; savemat(fp_zxyt_mat, 'd', "zc", adj_nu, 1, 0, (char *) temp_d_array, (char *) NULL); savemat(fp_zxyt_mat, 't', "z_label", 1, 1, 0, vgrid.label, (char *) NULL); savemat(fp_zxyt_mat, 't', "regrid_var", 1, 1, 0, get_name(regrid_variables, regrid_var), (char *) NULL); savemat(fp_zxyt_mat, 't', "regrid_method", 1, 1, 0, get_name(regrid_methods, regrid_method), (char *) NULL); savemat(fp_zxyt_mat, 't', "dbname", 1, 1, 0, DBNAME, (char *) NULL); savemat(fp_zxyt_mat, 't', "adcpsect_time", 1, 1, 0, time_stamp, (char *) NULL); temp_d_array[0] = (double) year_base; savemat(fp_zxyt_mat, 'd', "year_base", 1, 1, 0, (char *) temp_d_array, (char *) NULL); mat_zxyt_ptr = start_mat_array(fp_zxyt_mat, 'd', "xyt", 'c', 3); mat_uv_ptr = start_mat_array(fp_uv_mat, 'd', "uv", 'c', adj_nu); /* * ---> Start the loop through time ranges listed in the control file. */ print_header(fp_out); if (combined) initialize_ensemble(); while (get_time_range(fp_cnt, &start, &end_time) == 1) { check_dbsrch(TIME_SEARCH, (char *) &start); if (use_ctd) { check_dbset(CTD_DB, "DBSET, beginning of ctd block"); check_dbsrch(TIME_SEARCH, (char *) &start); get_ctd_profile(&ctd); printf("\n ctd position: %f %f", ctd.position.lat, ctd.position.lon); check_dbset(ADCP_DB, "DBSET, after get_ctd_profile"); } if (combined && !contour) print_time_range(fp_out, &start, &end_time); print_time_range(stdout, &start, &end_time); if (separate) initialize_ensemble(); /* * ---> Start the loop through profiles within the time range. */ while (!IERR && check_time(&adcp.time, &start, &end_time)) { if (verbose) { printf("\n"); write_ymdhms_time(stdout, &adcp.time); } adcp.ddtime = year_day(&adcp.time, year_base); if (single_ens) initialize_ensemble(); increment_ensemble(); DBMOVE(&STEPS, &IERR); /* IERR is checked at the start of the loop */ if (single_ens) end_ensemble(); } /* end of loop through profiles within a time * range */ if (separate) end_ensemble(); } /* end of loop through time ranges read from * control file */ puts("\n End of control file\n"); if (combined) end_ensemble(); time(&clock_end); run_time = (int) (clock_end - clock_start); fprintf(fp_out, "\n Elapsed time: %d:%02d (min:sec)\n", run_time / 60, run_time % 60); fclose(fp_out); if (contour) close_contour(); if (vector) fclose(fp_vec); if (transport) close_transport(); if (seq_cor) close_seqcor(); if (nav.method) close_navigation(); DBCLOSE(&IERR); check_error(IERR, "DBCLOSE"); if (use_ctd) { check_dbset(CTD_DB, "DBSET, before closing ctd"); DBCLOSE(&IERR); check_error(IERR, "DBCLOSE ctd"); } free_unistat(&u_stat); free_unistat(&v_stat); free_unistat(&adcp_pos_stat); free_unistat(&adcp_ship_stat); free_unistat(&tr_temp_stat); close_vgrid(); free(temp_d_array); end_mat_array(mat_zxyt_ptr); end_mat_array(mat_uv_ptr); fclose(fp_zxyt_mat); fclose(fp_uv_mat); printf("\nADCPSECT completed.\n"); return; }