/*************************************************************** FLUX.C module for ADCPSECT This module integrates in the horizontal, from ensemble to ensemble, using trapezoidal integration. If the input has already been integrated in the vertical by selecting the "integrate" regridding option, then the result is the transport. Tue 01-23-90 Changed time variable name from 'end' to (JR) 'end_time' because of problems on the Sun Tue 03-05-91 Added a transport direction option "along" which results in the line integral along the ship track being calculated: udx+vdy. ***************************************************************/ #include #include #include #include "dbext.h" #include "ioserv.h" #include "vector.h" #include "use_db.h" #include "vstat.h" #include "pos_to_m.h" #include "adcp.h" #include "adcpsect.h" #define DIR_ALONG 2 #define DIR_ALL 0 extern YMDHMS_TIME_TYPE start, end_time; /* for time range */ int transport = FALSE; /* global flag */ extern char out_filename[]; extern char time_stamp[]; static char tran_filename[NPATHMAX]; static FILE *fp_tran; static FILE *fp_tran_mat; static MAT_ARRAY_FP_TYPE *mat_ptr; /*--> for transport calculation using trapezoidal rule */ static double u_prev[NDMAX], v_prev[NDMAX], tran_inc[NDMAX]; static double u_this[NDMAX], v_this[NDMAX]; static double u_nan[NDMAX]; /* temporary array with BADFLOAT converted to NaN */ static DOUBLE_LL_TYPE pos_prev; static UNISTAT_TYPE u_tran, v_tran; static double tran_units = 1.0; static int tran_first_ensemble = 1; /* flag to suppress integration */ static int tran_min_npts; /* minimum number for transport output */ /* This threshold will determine whether integration proceeds at a given bin. If it stops, the sum is not zeroed or changed to bad; rather, the number of points in the integral is reduced. Hence, the user must look at N in the printout to determine whether all profiles were included in the integration at a given depth. */ static int tran_dir; /* direction of transport rel to track */ NAME_LIST_ENTRY_TYPE tran_directions[] = { {"right" , 1}, {"left" , -1}, {"all" , DIR_ALL}, {"along" , DIR_ALONG}, {NULL , 0} }; NAME_LIST_TYPE transport_name_lists[] = { tran_directions }; PARAMETER_LIST_ENTRY_TYPE transport_params[] = { {" minimum_npts=" , " %d" , &tran_min_npts , TYPE_INT}, {" units=" , " %lf" , &tran_units , TYPE_DOUBLE}, {" direction=" , " %s" , &tran_dir , TYPE_TRANS}, {NULL , NULL , NULL , 0} }; #if PROTOTYPE_ALLOWED int init_transport(FILE *fp_cnt) #else int init_transport(fp_cnt) FILE *fp_cnt; #endif { check_error(get_parameters(fp_cnt, transport_params, transport_name_lists), "transport parameters"); transport = TRUE; /* Make a transport output file name, and open the file. */ new_extension(tran_filename, out_filename, ".trn"); fp_tran = check_fopen(tran_filename, "w"); new_extension(tran_filename, out_filename, "_tr.mat"); fp_tran_mat = check_fopen(tran_filename, "wb"); savemat(fp_tran_mat, 't', "direction", 1, 1, 0, get_name(tran_directions, tran_dir), (char *)NULL); savemat(fp_tran_mat, 't', "adcpsect_flux_time", 1, 1, 0, time_stamp, (char *)NULL); return(0); } /* init_transport */ /* The allocation of arrays must be done in a separate function because it must wait until adj_ndepth is known. */ #if PROTOTYPE_ALLOWED void allocate_transport(int adj_nu) #else void allocate_transport(adj_nu) int adj_nu; #endif { allocate_unistat(&u_tran, adj_nu, "U transport", U_MEAN); allocate_unistat(&v_tran, adj_nu, "V transport", U_MEAN); /* Transport calculation will be cumulative for the whole run of adcpsect, so this is the only place that zeroing is needed. If malloc zeros the new memory, it would not really be needed even here. */ zero_unistat(&u_tran); zero_unistat(&v_tran); mat_ptr = start_mat_array(fp_tran_mat, 'd', "transport", 'c', adj_nu); } /* allocate_transport */ #if PROTOTYPE_ALLOWED void integrate_ds(double *last_ptr, double *now_ptr, int adj_nu, double distance, double *inc_ptr, int sign) #else void integrate_ds(last_ptr, now_ptr, adj_nu, distance, inc_ptr, sign) double *last_ptr, *now_ptr; int adj_nu, sign; double distance, *inc_ptr; #endif { double d, sum; int i; int last_pos, now_pos; /* On entry, sign is the desired sign of the product, u*ds. By reversing sign if ds < 0, we make sign into the desired sign of u alone. */ if (distance < 0.0) sign *= -1; distance *= 0.5; /* for trapezoidal rule */ for (i=0; i 0.0); now_pos = (*now_ptr > 0.0); sum = fabs(*last_ptr) + fabs(*now_ptr); if ( (sign > 0 && !last_pos && !now_pos) || (sign < 0 && last_pos && now_pos) ) *inc_ptr = 0.0; /* not in an integration region */ else if ( last_pos && !now_pos) { /* transition plus to minus */ if (sign > 0) /* ending an integration region */ { d = distance * fabs(*last_ptr) / sum; *inc_ptr = d * (*last_ptr); } else /* starting an integration region */ { d = distance * fabs(*now_ptr) / sum; *inc_ptr = d * (*now_ptr); } } else if ( !last_pos && now_pos) { /* transition minus to plus */ if (sign < 0) /* ending an integration region */ { d = distance * fabs(*last_ptr) / sum; *inc_ptr = d * (*last_ptr); } else /* starting an integration region */ { d = distance * fabs(*now_ptr) / sum; *inc_ptr = d * (*now_ptr); } } else /* entirely in an integration region */ { d = distance; *inc_ptr = d * (*last_ptr + *now_ptr); } } else /* sign == 0 */ { d = distance; *inc_ptr = d * (*last_ptr + *now_ptr); } } else *inc_ptr = BADFLOAT; } } /* integrate_ds() */ #if PROTOTYPE_ALLOWED void calculate_transport(UNISTAT_TYPE *u_stat, UNISTAT_TYPE *v_stat, UNISTAT_TYPE *adcp_pos_stat) #else void calculate_transport(u_stat, v_stat, adcp_pos_stat) UNISTAT_TYPE *u_stat, *v_stat, *adcp_pos_stat; #endif { double dx, dy; int i; memcpy(u_this, u_stat->sum, u_stat->n * sizeof(double)); memcpy(v_this, v_stat->sum, v_stat->n * sizeof(double)); if (!tran_first_ensemble) { dx = dlon_to_meters(adcp_pos_stat->sum[LONGITUDE] - pos_prev.lon, pos_prev.lat); dy = dlat_to_meters(adcp_pos_stat->sum[LATITUDE] - pos_prev.lat, pos_prev.lat); /* To block the integration if there are not enough points in the present ensemble, mark u and v as bad. */ for (i=0; in; i++) { if (u_stat->npts[i] < tran_min_npts) { u_this[i] = BADFLOAT; v_this[i] = BADFLOAT; } } if (tran_dir == DIR_ALONG) integrate_ds(u_prev, u_this, u_stat->n, dx, tran_inc, DIR_ALL); else integrate_ds(u_prev, u_this, u_stat->n, dy, tran_inc, tran_dir); update_unistat_d(&u_tran, tran_inc, u_stat->n); /* Since we are passing -dx, not dx, into integrate_ds for the v component, we pass tran_dir, not -tran_dir. */ if (tran_dir == DIR_ALONG) integrate_ds(v_prev, v_this, u_stat->n, dy, tran_inc, DIR_ALL); else integrate_ds(v_prev, v_this, u_stat->n, -dx, tran_inc, tran_dir); update_unistat_d(&v_tran, tran_inc, u_stat->n); } pos_prev.lon = adcp_pos_stat->sum[LONGITUDE]; pos_prev.lat = adcp_pos_stat->sum[LATITUDE]; memcpy(u_prev, u_this, u_stat->n * sizeof(double)); memcpy(v_prev, v_this, u_stat->n * sizeof(double)); tran_first_ensemble = 0; } /* calculate_transport() */ /*************************************************************** write_transport() (derived from write_unistat) The input variable d is a pointer to a float array used for labelling the rows of output. If the statistics are calculated from velocity profiles, then d is the array of depths. If fourier transformed profiles were used, then d is the array of vertical wavenumbers. If d is a NULL pointer, then the array index, starting from 1, will be printed. ***************************************************************/ #define depth (d==NULL ? (float) i+2 : d[i+1]) #if PROTOTYPE_ALLOWED void write_transport(FILE *fpt, UNISTAT_TYPE *u, UNISTAT_TYPE *v, double scale, float *d, char *d_name, int nd_prec) #else void write_transport(fpt, u, v, scale, d, d_name, nd_prec) FILE *fpt; UNISTAT_TYPE *u, *v; double scale; /* e.g., set scale to 0.01 to convert m/s to cm/s */ float *d; /* if it is a null pointer, then the index will be used */ char *d_name; /* for labelling the "depth" column */ int nd_prec; /* precision (decimal places) for "depth" column */ #endif { int i, np; double uu, vv; double uus = 0.0, vvs = 0.0; check_error( (scale == 0.0), "write_transport; scale is zero"); fprintf(fpt, "\nADCPSECT %s", time_stamp); fprintf(fpt, "\nTransports scaled by %e", scale); fprintf(fpt, "\nDepth, bin, or density variable gives lower boundary of integration"); fprintf(fpt, "\n --------horizontal--------"); fprintf(fpt, " --horizontal and vertical--"); if (tran_dir == DIR_ALONG) fprintf(fpt, "\n%12s N udx vdy sum udx vdy sum ", d_name); else fprintf(fpt, "\n%12s N udy -vdx sum udy -vdx sum ", d_name); for (i=0; in; i++) if (depth <= ADJ_BADFLOAT) { if ( (np = (u->npts)[i]) >0 ) { uu = (u->sum)[i]/scale; vv = (v->sum)[i]/scale; uus += uu; vvs += vv; fprintf(fpt, "\n%10.*f %6d %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f", nd_prec, depth, np, uu, vv, uu+vv, uus, vvs, uus+vvs); } } fprintf(fpt, "\n"); } /* write_transport() */ #undef depth #if PROTOTYPE_ALLOWED void write_tran_mat(MAT_ARRAY_FP_TYPE *mat_ptr, double *u, double *v, int n) #else void write_tran_mat(mat_ptr, u, v, n) MAT_ARRAY_FP_TYPE *mat_ptr; double *u, *v; int n; #endif { tonan_d(u_nan, u, n); add_to_mat_array(mat_ptr, (char *)u_nan); tonan_d(u_nan, v, n); add_to_mat_array(mat_ptr, (char *)u_nan); } #if PROTOTYPE_ALLOWED void print_transport(UNISTAT_TYPE *adcp_pos_stat, float *y, char *label, int n_zdigits) #else void print_transport(adcp_pos_stat, y, label, n_zdigits) UNISTAT_TYPE *adcp_pos_stat; float *y; char *label; int n_zdigits; #endif { extern YMDHMS_TIME_TYPE start, end_time; write_vgrid_params(fp_tran); print_time_range(fp_tran, &start, &end_time); fprintf(fp_tran, "\n ADCP mean position: %5.2f lat, %6.2f long\n", adcp_pos_stat->sum[LATITUDE], adcp_pos_stat->sum[LONGITUDE]); write_transport(fp_tran, &u_tran, &v_tran, tran_units, y, label, n_zdigits); write_tran_mat(mat_ptr, u_tran.sum, v_tran.sum, u_tran.n); fprintf(fp_tran, "\n\n"); } #if PROTOTYPE_ALLOWED void close_transport(void) #else void close_transport() #endif { fclose(fp_tran); end_mat_array(mat_ptr); fclose(fp_tran_mat); free_unistat(&u_tran); free_unistat(&v_tran); }