/*************************************************************** CTD module for ADCPSECT */ #include "geninc.h" #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 "adcpsect.h" /* includes ctd.h */ #include "data_id.h" #define TRUE 1 #define FALSE 0 #define sqr(x) (x)*(x) #define good_float(x) ( (x) < ADJ_BADFLOAT ) int use_ctd = FALSE; /* global flag; ctd data in use */ static FILE_NAME_TYPE ctd_dbname; /* ctd option */ double nominal_svel = NOMINAL_SVEL; #if PROTOTYPE_ALLOWED int init_ctd(FILE *fp_cnt) #else int init_ctd(fp_cnt) FILE *fp_cnt; #endif { use_ctd = TRUE; check_error( (fscanf(fp_cnt," dbname: %79s", ctd_dbname) != 1), "Can't read ctd_dbname"); check_dbopen(CTD_DB, ctd_dbname, READ_ONLY, DIR_IN_MEMORY); printf("\n\n DATABASE: %s SUCCESSFULLY OPENED",ctd_dbname); check_dbset(ADCP_DB, "DBSET, after opening ctd db"); return(0); } #if PROTOTYPE_ALLOWED void write_ctd_params(FILE *fp) #else void write_ctd_params(fp) FILE *fp; #endif { fprintf(fp, "\n CTD Database: %s", ctd_dbname); fprintf(fp, "\n depth correction with nominal sound speed = %7.2f", nominal_svel); } #if PROTOTYPE_ALLOWED void set_ctd_grid(CTD_PROFILE_TYPE *ctd, int regrid_var) #else void set_ctd_grid(ctd, regrid_var) CTD_PROFILE_TYPE *ctd; int regrid_var; #endif { if (regrid_var == RV_SVAN) ctd->density_var = CTD_SVAN; else if (regrid_var == RV_POT_SVAN) ctd->density_var = CTD_POT_SVAN; else if (regrid_var == RV_SIGMA_THETA) ctd->density_var = CTD_SIGMA_THETA; else if (regrid_var == RV_SIGMA_T) ctd->density_var = CTD_SIGMA_T; else ctd->density_var = CTD_NONE; } #if PROTOTYPE_ALLOWED void get_ctd_profile(CTD_PROFILE_TYPE *prof) #else void get_ctd_profile(prof) CTD_PROFILE_TYPE *prof; #endif { double sigma; /* dummy used in SVAN */ double theta; int step = 1; int ngaps = 0; /* max number of interpolations done */ int ierr; int n, ng; int i; get_latlon(&(prof->position)); check_dbget(TIME, (char *) &(prof->time), sizeof(prof->time), "DBGET, time"); prof->n = check_dbget_f(P, prof->pressure, NPMAX, "DBGET_F, pressure"); n = check_dbget_f(TEMPERATURE, prof->temperature, NPMAX, "DBGET_F, temperature"); prof->n = min_val(prof->n, n); n = check_dbget_f(SALINITY, prof->salinity, NPMAX, "DBGET_F, salinity"); prof->n = min_val(prof->n, n); /* It is assumed that the CTD data are nearly perfect, so this gap filling will only occasionally be called into play, and only for small gaps. */ ng = fill_gaps(prof->salinity, prof->n); ngaps = ng; ng = fill_gaps(prof->temperature, prof->n); ngaps = max_val(ngaps, ng); ng = fill_gaps(prof->pressure, prof->n); ngaps = max_val(ngaps, ng); if (ngaps == prof->n) prof->n = 0; /* don't try to use this profile */ for (i=0; i < prof->n; i++) { prof->depth[i] = p_depth( prof->pressure[i], prof->position.lat); prof->svel[i] = SVEL( prof->salinity[i], prof->temperature[i], prof->pressure[i]); switch (prof->density_var) { case CTD_SVAN: prof->sigma_sva[i] = SVAN( prof->salinity[i], prof->temperature[i], prof->pressure[i], &sigma); break; case CTD_SIGMA_THETA: theta = THETA(prof->salinity[i], prof->temperature[i], 0.0, prof->pressure[i]); SVAN( prof->salinity[i], theta, 0.0, &sigma); prof->sigma_sva[i] = sigma; break; case CTD_POT_SVAN: theta = THETA(prof->salinity[i], prof->temperature[i], 0.0, prof->pressure[i]); prof->sigma_sva[i] = SVAN( prof->salinity[i], theta, 0.0, &sigma); break; case CTD_SIGMA_T: SVAN( prof->salinity[i], prof->temperature[i], 0.0, &sigma); prof->sigma_sva[i] = sigma; break; case CTD_NONE: prof->sigma_sva[i] = BADFLOAT; break; default: printf("\n%d is an invalid density variable specification\n", prof->density_var); exit(-1); } } if (prof->n == 0) /* try the next one */ { printf("\nBad CTD profile."); write_ymdhms_time(stdout, &(prof->time)); printf(" position: %f %f\n", prof->position.lat, prof->position.lon); DBMOVE(&step, &ierr); check_error(ierr, "DBMOVE in get_ctd_profile"); get_ctd_profile(prof); /* recursive call */ } return; /* only exit */ } /* get_ctd_profile */ /*************************************************************** correct_depths() The ADCP depth correction factor, to be multiplied times the nominal depth, is the average sound speed from the transducer to the depth bin, divided by the nominal value used in the instrument. Averaging from the center of the first bin, as done here, is close enough. (After 89/02/25 the first depth bin is at zero meters, rather than the actual position of the first recorded ADCP bin. This is generally closer to the transducer depth than the center of the first ADCP bin is.) Input: a pointer to a ctd profile with sound speed, the array of depths to be corrected, and a pointer to the number of elements in the depth array. Output: the depth array is corrected in place, and the number of elements may be reduced if the ctd profile did not extend deep enough. With each call, this routine checks to see if the ctd profile is the same as the last one (same time). Only if it is different does the routine calculate the correction factor. */ #if PROTOTYPE_ALLOWED void correct_depths(CTD_PROFILE_TYPE *ctd, float *depth, int *nd_ptr) #else void correct_depths(ctd, depth, nd_ptr) CTD_PROFILE_TYPE *ctd; float *depth; int *nd_ptr; #endif { int n0, nout; int i; double svel0; static YMDHMS_TIME_TYPE last_time; static float depth_correction[NDMAX]; static int nd_ctd; /* 94/02/10 jr+ */ if (TIMDIF(&last_time, &(ctd->time))) /* true -> new ctd profile */ { LONG long_0 = 0L; DIFTIM(&(ctd->time), &last_time, &long_0); /* last equals present */ /* Calculate soundspeed at the nominal ADCP depths. */ regrid(ctd->depth, ctd->svel, depth, depth_correction, ctd->n, *nd_ptr, &n0, &nout); /* If many points are lost, something is probably very wrong. */ if( (n0 > 5) || ((n0 + nout) < (*nd_ptr-5))) { printf("\nWARNING: available depth range in correct_depths is"); printf("\n %f to %f, at ctd time ", depth[n0], depth[n0 + nout - 1] ); write_ymdhms_time(stdout, &(ctd->time)); } /* If there is a gap at the start, assume constant soundspeed */ for (i=0; i