/***************************************************************** FIX_TEMP.C Within a time range, this program corrects the transducer temperature (ancil1.tr_temp), which recorded incorrectly earler and those variables calculated from the transducer temperature including sound speed (ancil1.snd_spd_used), U, V, W and Error velocity. Set offset_temperature= 0 if input true_temperature will be used. Otherwise true_temperature = original_temp + offset_temperature will be applied in velocity correction. Where original_temp is the temperature original recorded in the ping data file. Willa Zhu 03-07-1990 Modified Mon 07-09-1990 by EF; new parameter original_soundspeed is the soundspeed that was input to the DAS if the DAS was not set to calculate soundspeed from temperature, or 0 otherwise. In the former case one can compensate for errors in manual input of soundspeed to the DAS; in the latter one compensates for errors in the soundspeed due to temperature errors. ------ Modified to add offset_temperature input parameter (WZ, 06/10/91) 92/01/15 (JR) Modified to allow user to provide temperature & salinity values from an external file with 2 columns: ADCP time (in decimal days) and the true temperature, a temperature offset, or salinity. 93/07/27 (JR) Moved n_vel initialization inside the while loop so all bins get grabbed and rotated even in the case where the conf.num_bins may change within the database 93/10/29 (JR) Added rescaling of bottom track velocities also. Previously, only U, V, W, and E were rescaled. 2000/06/23 (EF) Added ability to specify true_soundspeed= as a file. 2000/07/27 (EF) Changed MIN_SAL to 0, MAX_SAL to 40. ------------------------------------------------------------------ control file structure: dbname: (e.g. WEP) original_soundspeed= offset_temperature= true_temperature= true_salinity= year_base= end time_ranges: (list of YMDHMS time pairs:) (yy/mm/dd hh:mm:ss) to (yy/mm/dd hh:mm:ss) . . ****************************************************************/ #include "geninc.h" #include "ioserv.h" #include "use_db.h" #include "adcp.h" #include "data_id.h" #include "dpmask.h" #define sqr(x) (x)*(x) #define equal(x,y) (fabs((x)-(y)) < 0.0001) #define good_float(x) ( (x) < ADJ_BADFLOAT ) #define MAX_NBINS 128 #define OFFSET 1 #define TRUE_TEMP 2 #define TRUE_SALINITY 3 #define TRUE_SOUNDSPEED 4 #define MIN_TEMP -5.0 /* deg C */ #define MAX_TEMP 45.0 #define MIN_SAL 0.0 /* PSU */ #define MAX_SAL 40.0 #define MIN_SOUNDSPEED 1400 /* m/s */ #define MAX_SOUNDSPEED 1600 #define FV_UNUSED 0 #define FV_CONSTANT 1 #define FV_FILE 2 typedef struct { int source; double value; FILE *fp; } FILE_VALUE_TYPE; #if PROTOTYPE_ALLOWED int check_file(FILE *fp, RANGE_TYPE range, int year_base, int type); int read_param(FILE *fp_cnt, char *arg, int code); #endif #if PROTOTYPE_ALLOWED int read_param(FILE *fp_cnt, char *arg, int code) #else int read_param(fp_cnt, arg, code) FILE *fp_cnt; char *arg; int code; #endif { FILE_VALUE_TYPE *fv; FILE_NAME_TYPE filename; double value; fv = (FILE_VALUE_TYPE *) arg; if (fscanf(fp_cnt, "%s", filename) != 1) { fprintf(stderr, " ERROR: Unable to read input parameter from control file\n\n"); exit(-1); } if (sscanf(filename, " %lf", &value) != 1) /* user gave filename, not double value */ { if ( (fv->fp = fopen(filename, "rt")) == NULL ) { fprintf(stderr, " ERROR: Unable to open input file %s\n\n", filename); exit(-1); } if (ECHO) fprintf(stdout, " file: %s", filename); fv->source = FV_FILE; } else { fv->value = value; fv->source = FV_CONSTANT; if (ECHO) fprintf(stdout, " %f", fv->value); switch (code) /* check for reasonable values; in the case where values are supplied in a file, the check is done later together with the check for 1-1 match with ADCP times--see check_file() */ { case TRUE_TEMP: if (fv->value != 0.0) /* 0.0 allowed for backward compatibility */ { if (fv->value < MIN_TEMP || fv->value > MAX_TEMP) { fprintf(stderr, " ERROR: Unreasonable value for temperature %f\n\n", fv->value); exit(-1); } } break; case TRUE_SALINITY: if (fv->value < MIN_SAL || fv->value > MAX_SAL) { fprintf(stderr, " ERROR: Unreasonable value for salinity %f\n\n", fv->value); exit(-1); } break; case OFFSET: /* too much hassle checking for reasonable offset: would need to actually step thru every tr_temp in the database and see if result is reasonable */ break; case TRUE_SOUNDSPEED: break; default: break; } } return(0); } /* Main routine */ #if PROTOTYPE_ALLOWED void do_it(FILE *fp_cnt) #else void do_it(fp_cnt) FILE *fp_cnt; #endif { int i, IERR, STEP=1; unsigned int nvalues, n_vel, n_bt, ancil_size, n_time = sizeof(YMDHMS_TIME_TYPE); YMDHMS_TIME_TYPE prftime; YMDHMS_TIME_TYPE start, end, /* for time range */ this_time; /* current profile time */ ANCILLARY_1_TYPE ancil1; SHORT u[MAX_NBINS], v[MAX_NBINS], w[MAX_NBINS], e[MAX_NBINS]; BOTTOM_TRACK_TYPE bt; int const_snd_spd; /* flag: 1 if sound speed was specified to the DAS, 0 if it was calculated by the DAS. */ double old_snd_spd; double true_temp, corr_snd_spd, scale_factor; long int range_ofs; RANGE_TYPE range; char buff[80]; double filetime, dbtime; int is_offset; static FILE_NAME_TYPE dbname; static FILE_VALUE_TYPE temp, offset, sal, sndspd; static double orig_snd_spd; static int year_base; static OPTION_TYPE options[] = { { "end" , 0 , NULL , (char *) NULL}, { "dbname:" , TYPE_STRING , op_get_param , dbname }, { "year_base=" , TYPE_INT , op_get_param , (char *) &year_base}, { "original_soundspeed=", TYPE_DOUBLE , op_get_param , (char *) &orig_snd_spd}, { "offset_temperature=" , OFFSET , read_param , (char *) &offset}, { "true_temperature=" , TRUE_TEMP , read_param , (char *) &temp}, { "true_salinity=" , TRUE_SALINITY, read_param , (char *) &sal}, { "true_soundspeed=" , TRUE_SOUNDSPEED, read_param , (char *) &sndspd}, { NULL , 0 , NULL , (char *) NULL} }; /* options[] */ /***************************************************************/ sin(0.0); /* to force linking of float routines */ /* defaults: */ orig_snd_spd = 0; year_base = -99; dbname[0] = '\0'; is_offset = 0; /* true temperature or offset? */ temp.fp = sal.fp = offset.fp = sndspd.fp = NULL; temp.value = sal.value = offset.value = sndspd.value = 0.0; temp.source = sal.source = offset.source = sndspd.source = FV_UNUSED; fprintf(stdout, "\n Pass 1: Checking control file...\n"); if (execute_options(fp_cnt, options, ECHO) <= 0) exit(-1); #if 0 if (!( (temp.fp && temp.value == 0.0 && offset.fp == NULL && offset.value == 0.0) || (offset.fp && offset.value == 0.0 && temp.fp == NULL && temp.value == 0.0) || (temp.fp == NULL && temp.value != 0.0 && offset.fp == NULL && offset.value == 0.0) || (offset.fp == NULL && offset.value != 0.0 && temp.fp == NULL && temp.value == 0.0) ) || (sal.fp == NULL && sal.value == 0.0) ) #else if ( (sndspd.source && (temp.source || sal.source || offset.source)) || (sndspd.source == UNUSED && (sal.source == UNUSED || (temp.source && offset.source) || (!temp.source && !offset.source)))) #endif { fprintf(stderr, " ERROR: You must supply either true_temperature OR offset_temperature\n"); fprintf(stderr, " AND true_salinity; OR true_soundspeed alone\n\n"); exit(-1); } is_offset = (offset.fp || offset.value != 0.0); if (orig_snd_spd < 1400.0) const_snd_spd = 0; else const_snd_spd = 1; /* check_dbopen(1, dbname, READ_WRITE, DIR_IN_MEMORY); */ check_dbopen(1, dbname, READ_WRITE, DIR_ON_DISK); if (year_base == -99) /* use first profile time as default */ { if (dbget_cnf(TIME, (char *) &prftime, &n_time, "getting time")) exit(-1); year_base = prftime.year; } check_error( (fscanf(fp_cnt, " time_ranges:") != 0), "reading time range start"); range.type = TIME_RANGE; range_ofs = ftell(fp_cnt); /* note down for processing pass later */ while (get_time_range(fp_cnt, &(range.ru.time.start), &(range.ru.time.end)) == 1) { switch(check_range(&range)) { case GAP: print_range(stdout, &range, "\n WARNING: Gap before range "); break; case OVERLAP: print_range(stdout, &range, "\n ERROR: Range overlaps preceding => "); exit(-1); case BAD_RANGE: print_range(stdout, &range, "\n ERROR: Bad range => "); exit(-1); } if (sndspd.fp && check_file(sndspd.fp, range, year_base, TRUE_SOUNDSPEED)) exit(-1); if (sal.fp && check_file(sal.fp, range, year_base, TRUE_SALINITY)) exit(-1); if (temp.fp && check_file(temp.fp, range, year_base, TRUE_TEMP)) exit(-1); else if (offset.fp && check_file(offset.fp, range, year_base, OFFSET)) exit(-1); } if (check_for_gap_at_end(&range) == GAP) { print_range(stdout, &range, "\n WARNING: Gap after range "); } /* rewind for pass 2: */ fseek(fp_cnt, range_ofs, SEEK_SET); if (sal.fp) fseek(sal.fp, 0L, SEEK_SET); if (temp.fp) fseek(temp.fp, 0L, SEEK_SET); else if (offset.fp) fseek(offset.fp, 0L, SEEK_SET); if (sndspd.fp) fseek(sndspd.fp, 0L, SEEK_SET); fprintf(stdout, "\n Pass 2: Updating database...\n"); /* ---> Start the loop through time ranges listed in the control file. */ while (get_time_range(fp_cnt, &start, &end) == 1) { check_dbsrch(TIME_SEARCH, (char *)&start); print_time_range(stdout, &start, &end); IERR = 0; /* ---> Start the loop through profiles within the time range. */ ancil_size = sizeof(ancil1); while (!IERR && check_time(&this_time, &start, &end)) { dbtime = year_day(&this_time, year_base); n_vel = MAX_NBINS * sizeof(SHORT); n_bt = sizeof(BOTTOM_TRACK_TYPE); if (sndspd.source == FV_FILE) while (getline_nc(sndspd.fp, buff, 80) != EOF && sscanf(buff, " %lf %lf", &filetime, &sndspd.value) == 2 && !equal(filetime, dbtime)); if (sal.source == FV_FILE) while (getline_nc(sal.fp, buff, 80) != EOF && sscanf(buff, " %lf %lf", &filetime, &sal.value) == 2 && !equal(filetime, dbtime)); if (temp.source == FV_FILE) while (getline_nc(temp.fp, buff, 80) != EOF && sscanf(buff, " %lf %lf", &filetime, &temp.value) == 2 && !equal(filetime, dbtime)); else if (offset.source == FV_FILE) while (getline_nc(offset.fp, buff, 80) != EOF && sscanf(buff, " %lf %lf", &filetime, &offset.value) == 2 && !equal(filetime, dbtime)); if ((IERR = dbget_cnf(ANCILLARY_1, (char *)&ancil1, &ancil_size, "getting ANCILARY_1")) == 0 && (IERR = dbget_cnf(U, (char *)u, &n_vel,"getting U")) == 0 && (IERR = dbget_cnf(V, (char *)v, &n_vel,"getting V")) == 0 && (IERR = dbget_cnf(W, (char *)w, &n_vel,"getting W")) == 0 && (IERR = dbget_cnf(ERROR_VEL,(char *)e, &n_vel, "getting ERROR_VEL")) == 0 && (IERR = dbget_cnf(BOTTOM_TRACK, (char *)&bt, &n_bt,"getting BOTTOM TRACK")) == 0) { /* 2000/12/12 PJ * Added a message in case n_vel becomes zero after one of the * calls above. This is unfortunately the case sometimes. */ if(n_vel == 0){ char str[299]; sprintf(str, "\nBad number of elements in profile at time %12.6f\n", dbtime); report_msg(str); } if (offset.source) true_temp = ancil1.tr_temp + offset.value; else true_temp = temp.value; if (sndspd.source) { corr_snd_spd = sndspd.value; } else { corr_snd_spd = 1449.2 + 4.6 * true_temp - 0.055 * true_temp*true_temp + 0.00029 * (true_temp*true_temp*true_temp) + (1.34 - 0.01*true_temp) * (sal.value-35.0); } /* one more check for reasonable values: */ if (corr_snd_spd < MIN_SOUNDSPEED || corr_snd_spd > MAX_SOUNDSPEED) { fprintf(stderr, " ERROR: Unreasonable value (%f) for soundspeed for dbtime %f\n\n", corr_snd_spd, dbtime); goto close; } if (const_snd_spd) old_snd_spd = orig_snd_spd; else old_snd_spd = ancil1.snd_spd_used; nvalues = n_vel / sizeof(SHORT); scale_factor = corr_snd_spd/old_snd_spd; for (i = 0; i < nvalues; i++) { if (u[i] != BADSHORT) u[i] = (double)u[i] * scale_factor; if (v[i] != BADSHORT) v[i] = (double)v[i] * scale_factor; if (w[i] != BADSHORT) w[i] = (double)w[i] * scale_factor; if (e[i] != BADSHORT) e[i] = (double)e[i] * scale_factor; } if (bt.u < ADJ_BADFLOAT) { bt.u = (double) bt.u * scale_factor; bt.v = (double) bt.v * scale_factor; } ancil1.snd_spd_used = corr_snd_spd; ancil1.tr_temp = true_temp; ancil_size = sizeof(ancil1); if((IERR = dbput_cnf(ANCILLARY_1, (char *) &ancil1, &ancil_size, "putting ANCILLARY_1")) == 0) if ((IERR = dbput_cnf(U, (char *)u, &n_vel, "putting U velocity")) == 0) if ((IERR = dbput_cnf(V, (char *)v, &n_vel, "putting V velocity")) == 0) if ((IERR = dbput_cnf(W, (char *)w, &n_vel, "putting W velocity")) == 0) if((IERR = dbput_cnf(ERROR_VEL, (char *) e, &n_vel, "putting ERROR_VEL")) == 0) if((IERR = dbput_cnf(BOTTOM_TRACK, (char *) &bt, &n_bt, "putting BOTTOM_TRACK")) == 0) if ( (IERR = set_block_dpmask(SOUND_SPEED_VEL_COR)) == 0) DBMOVE(&STEP, &IERR); } } /* end of loop through profiles within a time range */ } /* end of loop through time ranges read from control file */ close: DBCLOSE(&IERR); check_error(IERR, "DBCLOSE"); printf("\n fix_temp completed.\n"); return; } #if PROTOTYPE_ALLOWED int check_file(FILE *fp, RANGE_TYPE range, int year_base, int type) #else int check_file(fp, range, year_base, type) FILE *fp; RANGE_TYPE range; int year_base; int type; /* TRUE_SOUNDSPEED, TRUE_TEMP, OFFSET or TRUE_SALINITY */ #endif { YMDHMS_TIME_TYPE current; char buff[80]; double dbtime, filetime, value; ANCILLARY_1_TYPE anc1; if (goto_start_of_range(&range)) return(1); while (in_range((char *) ¤t, &range)) { dbtime = year_day(¤t, year_base); while (getline_nc(fp, buff, 80) != EOF && sscanf(buff, " %lf %lf", &filetime, &value) == 2 && !equal(filetime, dbtime)); /* match times */ if (feof(fp)) { fprintf(stderr, " ERROR: No file time to match database time %f\n\n", dbtime); return(-1); } /* check for reasonable value */ switch (type) { case OFFSET: check_dbget(ANCILLARY_1, (char *) &anc1, sizeof(ANCILLARY_1_TYPE), "getting ANCILLARY_1"); value += anc1.tr_temp; case TRUE_TEMP: if (value < MIN_TEMP || value > MAX_TEMP) { fprintf(stderr, " ERROR: Unreasonable temperature for time %f in file\n\n", filetime); exit(-1); } break; case TRUE_SALINITY: if (value < MIN_SAL || value > MAX_SAL) { fprintf(stderr, " ERROR: Unreasonable salinity for time %f in file\n\n", filetime); exit(-1); } break; case TRUE_SOUNDSPEED: if (value < MIN_SOUNDSPEED || value > MAX_SOUNDSPEED) { fprintf(stderr, " ERROR: Unreasonable value (%f) for soundspeed at time %f in file\n\n", value, filetime); exit(-1); } break; default: break; } switch(dbmove_cnf(1)) { case 0: break; case SEARCH_BEYOND_END: return(0); default: return(-1); } } return(0); }