diff --git a/Wave_gpu.cu b/Wave_gpu.cu index 4d6acf3..d2ee648 100644 --- a/Wave_gpu.cu +++ b/Wave_gpu.cu @@ -1753,7 +1753,8 @@ int main(int argc, char **argv) if (XParam.wavebndtype == 3) { //Reuse XBeach_GPU boundary. Similar to XBeach but in netcdf format and all self explanatory - //Not implemented yet + XParam = read_reuse_bndnc_head(XParam); + wavebnd = read_reuse_bndnc_vec(XParam); } if (XParam.wavebndtype == 4) { diff --git a/Wavestep.cu b/Wavestep.cu index a6b972f..cf872db 100644 --- a/Wavestep.cu +++ b/Wavestep.cu @@ -57,7 +57,7 @@ XBGPUParam waveinitGPU(XBGPUParam Param, std::vector wavebnd) } if (Param.wavebndtype >= 2) { - nwavbnd = ceil(Param.rtlength / Param.dtbc); + nwavbnd = ceil(Param.rtlength / Param.dtbc)+1; // +1 needed here } theta = (DECNUM *)malloc(ntheta*sizeof(DECNUM)); @@ -126,7 +126,7 @@ XBGPUParam waveinitGPU(XBGPUParam Param, std::vector wavebnd) if (Param.wavebndtype == 3) { // Reuse XBeach_GPU style wave boundary. same as normal XBeach but as a self documented netcdf file - //Not implemented yet + read_reuse_bndnc(Param, 0, Trep, qfile, Stfile); } if (Param.wavebndtype == 4) { @@ -144,7 +144,61 @@ XBGPUParam waveinitGPU(XBGPUParam Param, std::vector wavebnd) //void GenWGnLBW(XBGPUParam Param, int nf, int ndir, double * HRfreq, double * HRdir, double * HRSpec, float Trep, double * qfile, double * Stfile) GenWGnLBW(Param, nfHR, ndHR, HRfreq, HRdir, HRSpec, Trep, qfile, Stfile); //create2dnc(nfHR, ndHR, HRfreq[1] - HRfreq[0], HRdir[1] - HRdir[0], 0.0, HRfreq, HRdir, HRSpec); + + ////////////////////////////////////// + //Save to Netcdf file + ////////////////////////////////////// + double * yyfx, *ttfx, *thetafx; + double * qxtemp, *qytemp, *eetemp; + int tslenbc = nwavbnd; + + qxtemp = (double *)malloc(ny*tslenbc*sizeof(double)); + qytemp = (double *)malloc(ny*tslenbc*sizeof(double)); + eetemp = (double *)malloc(ny*Param.ntheta*tslenbc*sizeof(double)); + + yyfx = (double *)malloc(ny*sizeof(double)); + ttfx = (double *)malloc(tslenbc*sizeof(double)); + thetafx = (double *)malloc(Param.ntheta*sizeof(double)); + + for (int j = 0; j < Param.ny; j++) + { + yyfx[j] = j*Param.dx; + } + + for (int m = 0; m < tslenbc; m++) + { + ttfx[m] = m*Param.dtbc; + } + + for (int itheta = 0; itheta < Param.ntheta; itheta++) + { + thetafx[itheta] = itheta*(Param.dtheta) + Param.thetamin + 0.5f*Param.dtheta; + } + + //for Debugging + //create2dnc(ny, tslen, Param.dx, dtin, 0.0, yyfx, tin, qx); + //create3dnc(ny, Param.ntheta, tslen, Param.dx, Param.dtheta, dtin, 0.0, yyfx, thetafx, tin, zeta); + + for (int j = 0; j < Param.ny; j++) + { + for (int m = 0; m < tslenbc; m++) + { + qxtemp[m + j*tslenbc] = qfile[j + 0 * ny + m*ny * 4]; + qytemp[m + j*tslenbc] = qfile[j + 1 * ny + m*ny * 4]; + + for (int itheta = 0; itheta < Param.ntheta; itheta++) + { + eetemp[m + j*tslenbc + itheta*ny*tslenbc] = Stfile[j + itheta*ny + m*ny*Param.ntheta]; + } + } + } + + + createbndnc(tslenbc, ny, Param.ntheta, Param.dx, Param.dtheta, 0.0, ttfx, yyfx, thetafx, eetemp, qxtemp, qytemp); //Trep = 15.0;// + free(HRSpec); + free(HRfreq); + free(HRdir); @@ -382,10 +436,90 @@ void wavebnd(XBGPUParam Param, std::vector wavebndvec) if (Param.wavebndtype == 2) { - //Read new STfile and qfile + //Read new STfile and qfile XBeach style readXbbndstep(Param, wavebndvec, WAVstepinbnd - 1, Trep, qfile, Stfile); } + + if (Param.wavebndtype == 3) + { + // Reuse XBeach_GPU style wave boundary. same as normal XBeach but as a self documented netcdf file + read_reuse_bndnc(Param, WAVstepinbnd - 1, Trep, qfile, Stfile); + } + + if (Param.wavebndtype == 4) + { + //JONSWAP + //First generate a Highres 2D spec + double * HRfreq; + double * HRdir; + double * HRSpec; + + int nfHR, ndHR; + + makjonswap(Param, wavebndvec, WAVstepinbnd - 1, nfHR, ndHR, HRfreq, HRdir, HRSpec); + + //Then generate wave group timeseries based on that spectra + //void GenWGnLBW(XBGPUParam Param, int nf, int ndir, double * HRfreq, double * HRdir, double * HRSpec, float Trep, double * qfile, double * Stfile) + GenWGnLBW(Param, nfHR, ndHR, HRfreq, HRdir, HRSpec, Trep, qfile, Stfile); + ////////////////////////////////////// + //Save to Netcdf file + ////////////////////////////////////// + + // Stfile is not ordered teh way we want to save it to file so we need a temporary storage to rearange + double * yyfx, *ttfx, *thetafx; + double * qxtemp, *qytemp, *eetemp; + int tslenbc = (int)ceil(Param.rtlength / Param.dtbc)+1; + + qxtemp = (double *)malloc(ny*tslenbc*sizeof(double)); + qytemp = (double *)malloc(ny*tslenbc*sizeof(double)); + eetemp = (double *)malloc(ny*Param.ntheta*tslenbc*sizeof(double)); + + yyfx = (double *)malloc(ny*sizeof(double)); + ttfx = (double *)malloc(tslenbc*sizeof(double)); + thetafx = (double *)malloc(Param.ntheta*sizeof(double)); + + for (int j = 0; j < Param.ny; j++) + { + yyfx[j] = j*Param.dx; + } + + for (int m = 0; m < tslenbc; m++) + { + ttfx[m] = m*Param.dtbc; + } + + for (int itheta = 0; itheta < Param.ntheta; itheta++) + { + thetafx[itheta] = itheta*(Param.dtheta) + Param.thetamin + 0.5f*Param.dtheta; + } + + //for Debugging + //create2dnc(ny, tslen, Param.dx, dtin, 0.0, yyfx, tin, qx); + //create3dnc(ny, Param.ntheta, tslen, Param.dx, Param.dtheta, dtin, 0.0, yyfx, thetafx, tin, zeta); + + for (int j = 0; j < Param.ny; j++) + { + for (int m = 0; m < tslenbc; m++) + { + qxtemp[m + j*tslenbc] = qfile[j + 0 * ny + m*ny * 4]; + qytemp[m + j*tslenbc] = qfile[j + 1 * ny + m*ny * 4]; + + for (int itheta = 0; itheta < Param.ntheta; itheta++) + { + eetemp[m + j*tslenbc + itheta*ny*tslenbc] = Stfile[j + itheta*ny + m*ny*Param.ntheta]; + } + } + } + + + writebndnc(tslenbc, ny, Param.ntheta, Param.dx, Param.dtheta, wavebndvec[WAVstepinbnd-1].time, ttfx, yyfx, thetafx, eetemp, qxtemp, qytemp); + //Trep = 15.0;// + free(HRSpec); + free(HRfreq); + free(HRdir); + + } } if (Param.wavebndtype == 1) diff --git a/XBeachGPU.h b/XBeachGPU.h index 9781258..a8e194b 100644 --- a/XBeachGPU.h +++ b/XBeachGPU.h @@ -217,6 +217,8 @@ extern "C" void read3Dnc(int nx, int ny,int ntheta,char ncfile[],DECNUM * &ee); extern "C" void read2Dnc(int nx, int ny,char ncfile[],DECNUM * &hh); void createbndnc(int tslen, int ny, int ntheta, double dy, double dtheta, double totaltime, double * timevec, double *yy, double *theta, double * ee, double * qx, double * qy); +void writebndnc(int tslen, int ny, int ntheta, double dy, double dtheta, double totaltime, double * timevec, double *yy, double *theta, double * ee, double * qx, double * qy); +void read_reuse_bndnc(XBGPUParam Param, int rec, float &Trep, double * &qfile, double * &Stfile); void GenCstWave(XBGPUParam Param, std::vector wavebnd, float * theta, double * &Stfile, double * &qfile, double * &Tpfile); @@ -229,6 +231,9 @@ extern "C" void readbndhead(const char * wavebndfile,DECNUM &thetamin,DECNUM &th std::vector ReadJSWPBnd(XBGPUParam XParam); +XBGPUParam read_reuse_bndnc_head(XBGPUParam Param); +std::vector read_reuse_bndnc_vec(XBGPUParam Param); + //Below is for the new CPU routine extern "C" int mminusC(int ix,int nx); extern "C" int pplusC(int ix, int nx); @@ -285,6 +290,7 @@ std::vector readWNDfile(std::string WNDfilename, double grdalpha); std::vector ReadCstBnd(XBGPUParam XParam); double interptime(double next, double prev, double timenext, double time); double interp1D(int nx, double *x, double *y, double xx); +double interp1DMono(int nx, double *x, double *y, double xx); double Interp2(int nx, int ny, double *x, double *y, double *z, double xx, double yy); extern "C" void readbathyHead(std::string filename, int &nx, int &ny, double &dx, double &grdalpha); diff --git a/makjonswap.cpp b/makjonswap.cpp index 0034aa1..31a0bc0 100644 --- a/makjonswap.cpp +++ b/makjonswap.cpp @@ -83,7 +83,8 @@ void makjonswap(XBGPUParam Param, std::vector wavebnd, int step, i double scoeff = wavebnd[step].s; - //printf("fp=%f\n",fp); + printf("Generating JONSWAP spectrum: Hs=%f, Tp=%f, Dp=%f, gam=%f, s=%f\n",Hs,Tp,Dp,gam,scoeff); + write_text_to_log_file("Generating JONSWAP spectrum: Hs=" + std::to_string(Hs) + " Tp=" + std::to_string(Tp) + " Dp=" + std::to_string(Dp) + " gam=" + std::to_string(gam) + " s=" + std::to_string(scoeff) + ";"); //// double fnyq = 3.0f*fp; @@ -222,6 +223,11 @@ void makjonswap(XBGPUParam Param, std::vector wavebnd, int step, i void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir, double * HRSpec, float &Trep, double * &qfile, double * &Stfile) { + // Generating Boundary condition: Energy from wave group and Long bound waves + + printf("Generating Boundary condition: Energy from wave group and Long bound waves.\n"); + write_text_to_log_file("Generating Boundary condition: Energy from wave group and Long bound waves"); + int ny = Param.ny; int K; //Number of wave components @@ -414,6 +420,8 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir //} //interp1D(double *x, double *y, double xx) //thetagen[i] = interptime(HRdir[dprev + 1], HRdir[dprev], dtheta, dtheta*((number - cdf[dprev]) / (cdf[dprev + 1] - cdf[dprev]))); + + //Cannot use interp1DMono here because cdf is not monotonic thetagen[i] = interp1D(ndir, cdf, HRdir, number); //printf("thetagen[i]=%f\n", thetagen[i]); @@ -470,7 +478,7 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir //First assume that internal and bc - writing time step is the same double dtin = Param.dtbc; - int tslenbc = ceil(Param.rtlength / Param.dtbc);// (int)(Param.rtlength / Param.dtbc) + 1; + int tslenbc = ceil(Param.rtlength / Param.dtbc) + 1;// (int)(Param.rtlength / Param.dtbc) + 1; //Check whether the internal frequency is high enough to describe the highest frequency //wave train returned from frange(which can be used in the boundary conditions) @@ -565,7 +573,7 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir double dummy; for (int i = 0; i < K; i++) { - dummy = interp1D(nf, HRfreq, Sf, fgen[i]); + dummy = interp1DMono(nf, HRfreq, Sf, fgen[i]); vargenq[i] = min(vargen[i] ,dummy); } @@ -944,7 +952,7 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir //interpolate to boundary timeseries for (int m = 0; m < tslenbc; m++) { - Stfile[j + itheta*ny + m*ny*Param.ntheta] = interp1D(tslen, tin, E_tdir, m*Param.dtbc); + Stfile[j + itheta*ny + m*ny*Param.ntheta] = interp1DMono(tslen, tin, E_tdir, m*Param.dtbc); } } @@ -1323,63 +1331,15 @@ void GenWGnLBW(XBGPUParam Param, int nf, int ndir,double * HRfreq,double * HRdir for (int m = 0; m < tslenbc; m++) { - qfile[j + 0 * ny + m*ny * 4] = interp1D(tslen, tin, qtempx, m*Param.dtbc); - qfile[j + 1 * ny + m*ny * 4] = interp1D(tslen, tin, qtempy, m*Param.dtbc); + qfile[j + 0 * ny + m*ny * 4] = interp1DMono(tslen, tin, qtempx, m*Param.dtbc); + qfile[j + 1 * ny + m*ny * 4] = interp1DMono(tslen, tin, qtempy, m*Param.dtbc); qfile[j + 2 * ny + m*ny * 4] = 0.0; qfile[j + 3 * ny + m*ny * 4] = 0.0;// qtot[j + m*ny]; } } - ////////////////////////////////////// - //Save to Netcdf file - ////////////////////////////////////// - double * yyfx, *ttfx, *thetafx; - double * qxtemp, *qytemp, *eetemp; - - qxtemp = (double *)malloc(ny*tslenbc*sizeof(double)); - qytemp = (double *)malloc(ny*tslenbc*sizeof(double)); - eetemp = (double *)malloc(ny*Param.ntheta*tslenbc*sizeof(double)); - - yyfx = (double *)malloc(ny*sizeof(double)); - ttfx = (double *)malloc(tslenbc*sizeof(double)); - thetafx = (double *)malloc(Param.ntheta*sizeof(double)); - for (int j = 0; j < Param.ny; j++) - { - yyfx[j] = j*Param.dx; - } - - for (int m = 0; m < tslenbc; m++) - { - ttfx[m] = m*Param.dtbc; - } - - for (int itheta = 0; itheta < Param.ntheta; itheta++) - { - thetafx[itheta] = itheta*(Param.dtheta) + Param.thetamin + 0.5f*Param.dtheta; - } - - //for Debugging - create2dnc(ny, tslen, Param.dx, dtin, 0.0, yyfx, tin, qx); - create3dnc(ny, Param.ntheta, tslen, Param.dx, Param.dtheta, dtin, 0.0, yyfx, thetafx, tin, zeta); - - for (int j = 0; j < Param.ny; j++) - { - for (int m = 0; m < tslenbc; m++) - { - qxtemp[m + j*tslenbc] = qfile[j + 0 * ny + m*ny * 4]; - qytemp[m + j*tslenbc] = qfile[j + 1 * ny + m*ny * 4]; - - for (int itheta = 0; itheta < Param.ntheta; itheta++) - { - eetemp[m + j*tslenbc + itheta*ny*tslenbc] = Stfile[j + itheta*ny + m*ny*Param.ntheta]; - } - } - } - - - createbndnc(tslenbc, ny, Param.ntheta, Param.dx, Param.dtheta, 0.0, ttfx, yyfx,thetafx, eetemp, qxtemp, qytemp); ////////////////////////////////////// //Clean up and desallocate all that memory diff --git a/read_input.cpp b/read_input.cpp index 5093a40..f498e0d 100644 --- a/read_input.cpp +++ b/read_input.cpp @@ -21,6 +21,7 @@ using DECNUM = float; + XBGPUParam readXbbndhead(XBGPUParam Param) { std::ifstream fs(Param.wavebndfile); @@ -107,6 +108,8 @@ std::vector readXbbndfile(XBGPUParam Param) return wavebndvec; } + + extern "C" void readbndhead(const char * wavebnd, DECNUM &thetamin, DECNUM &thetamax, DECNUM &dtheta, DECNUM &dtwavbnd, int &nwavbnd) { FILE * fwav; diff --git a/tools.cpp b/tools.cpp index 0b286f9..29cc687 100644 --- a/tools.cpp +++ b/tools.cpp @@ -27,6 +27,8 @@ double interptime(double next, double prev, double timenext, double time) double interp1D(int nx, double *x, double *y, double xx) { + // non monotonic 1D interpolation + // This can be pretty slow so use interp1DMono for faster interpolation if x is monotonic (i.e. the spacing is constant) double yy; double prevx = x[0]; double nextx = x[1]; @@ -58,6 +60,36 @@ double interp1D(int nx, double *x, double *y, double xx) yy = prevy + (xx - prevx) / (nextx - prevx)*(nexty - prevy); return yy; } + +double interp1DMono(int nx, double *x, double *y, double xx) +{ + // interpolation if x is monotonically increasing (i.e. the spacing is constant) + double yy; + double prevx = x[0]; + double nextx = x[1]; + double prevy = y[0]; + double nexty = y[1]; + double dx = nextx - prevx; + int indx = 0; + int indxp; + + double diffx = max(xx-x[0],0.0);//This has to be positive! + + indx = (int)floor(diffx / dx); + + + + indxp = (int)min(indx*1.0 + 1, nx*1.0 - 1); + prevx = x[indx]; + nextx = x[indxp]; + prevy = y[indx]; + nexty = y[indxp]; + + + yy = prevy + (xx - prevx) / (nextx - prevx)*(nexty - prevy); + return yy; +} + double Interp2(int nx, int ny, double *x, double *y, double *z, double xx, double yy) //double BilinearInterpolation(double q11, double q12, double q21, double q22, double x1, double x2, double y1, double y2, double x, double y) { diff --git a/writenetcdf.cpp b/writenetcdf.cpp index 7ef41f6..ca31726 100644 --- a/writenetcdf.cpp +++ b/writenetcdf.cpp @@ -22,6 +22,7 @@ using DECNUM = float; void handle_error(int status) { if (status != NC_NOERR) { fprintf(stderr, "Netcdf %s\n", nc_strerror(status)); + write_text_to_log_file("Netcdf " + std::string(nc_strerror(status))); //fprintf(logfile, "Netcdf: %s\n", nc_strerror(status)); exit(-1); } @@ -991,4 +992,255 @@ void createbndnc(int tslen, int ny, int ntheta, double dy, double dtheta, double +} + +void writebndnc(int tslen, int ny, int ntheta, double dy, double dtheta, double totaltime, double * timevec, double *yy, double *theta, double * ee, double * qx, double * qy) +{ + int status; + int ncid, xx_dim, yy_dim, time_dim, p_dim, ee_var_id, qx_var_id, qy_var_id, recid; + + size_t nxx, nyy, ntt; + static size_t nrec; + static size_t start_ee[] = { 0, 0, 0, 0 }; // start at first value + static size_t count_ee[] = { 1, ntheta, ny, tslen }; + static size_t start_q[] = { 0, 0, 0 }; // start at first value + static size_t count_q[] = { 1, ny, tslen }; + int time_id, xx_id, yy_id, tt_id; // + nxx = tslen; + nyy = ny; + ntt = ntheta; + + //create the netcdf dataset + status = nc_open("XBG_bnd_reuse.nc", NC_WRITE, &ncid); + + //read id from time dimension + status = nc_inq_unlimdim(ncid, &recid); + status = nc_inq_dimlen(ncid, recid, &nrec); + //printf("nrec=%d\n",nrec); + + //read file for variable ids + status = nc_inq_varid(ncid, "time", &time_id); + status = nc_inq_varid(ncid, "ee_bnd", &ee_var_id); + status = nc_inq_varid(ncid, "qx_bnd", &qx_var_id); + status = nc_inq_varid(ncid, "qy_bnd", &qy_var_id); + + start_ee[0] = nrec;// + start_q[0] = nrec;// + + + static size_t tst[] = { nrec }; + + + + //Provide values for variables + status = nc_put_var1_double(ncid, time_id, tst, &totaltime); + + + status = nc_put_vara_double(ncid, ee_var_id, start_ee, count_ee, ee); + status = nc_put_vara_double(ncid, qx_var_id, start_q, count_q, qx); + status = nc_put_vara_double(ncid, qy_var_id, start_q, count_q, qy); + + + //close file + status = nc_close(ncid); + + + +} + +XBGPUParam read_reuse_bndnc_head(XBGPUParam Param) +{ + //read the dimentions of grid, levels and time + int status; + int ncid, recid, theta_dimid, theta_varid, tt_dimid, tt_varid; + + static size_t nrec, ntheta, ntt; + status = nc_open(Param.wavebndfile.c_str(), NC_NOWRITE, &ncid); + + //read id from time dimension + status = nc_inq_unlimdim(ncid, &recid); + status = nc_inq_dimlen(ncid, recid, &nrec); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimid(ncid, "theta", &theta_dimid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimlen(ncid, theta_dimid, &ntheta); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_varid(ncid, "theta", &theta_varid); + if (status != NC_NOERR) handle_error(status); + + double thetamin, thetamax, dtheta; + + + size_t start[] = { 0 }; + size_t count[] = { 1 }; + status = nc_get_vara_double(ncid, theta_varid, start, count, &thetamin); + start[0] = ntheta-1; + status = nc_get_vara_double(ncid, theta_varid, start, count, &thetamax); + + Param.ntheta = ntheta; + Param.thetamin = thetamin; + Param.thetamax = thetamax; + + Param.dtheta = (Param.thetamax - Param.thetamin) / Param.ntheta; + + status = nc_inq_dimid(ncid, "tt", &tt_dimid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimlen(ncid, tt_dimid, &ntt); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_varid(ncid, "tt", &tt_varid); + if (status != NC_NOERR) handle_error(status); + + double rtl,rtlp; + + start[0] = ntt - 1; + status = nc_get_vara_double(ncid, tt_varid, start, count, &rtl); + start[0] = ntt - 2; + status = nc_get_vara_double(ncid, tt_varid, start, count, &rtlp); + + + //close file + status = nc_close(ncid); + + Param.rtlength = ntt*(rtl - rtlp); // Not totally sure why this is bigger than rtl... + Param.dtbc = rtl - rtlp; + + return Param; +} + +std::vector read_reuse_bndnc_vec(XBGPUParam Param) +{ + std::vector wavebndvec; + + Wavebndparam wavebndline; + + //read the dimentions of grid, levels and time + int status; + int ncid, recid, theta_dimid, theta_varid, tt_dimid, tt_varid; + + static size_t nrec, ntheta, ntt; + status = nc_open(Param.wavebndfile.c_str(), NC_NOWRITE, &ncid); + + //read id from time dimension + status = nc_inq_unlimdim(ncid, &recid); + status = nc_inq_dimlen(ncid, recid, &nrec); + if (status != NC_NOERR) handle_error(status); + + //close file + status = nc_close(ncid); + + //rtlength has just been calculated by teh previous function + + for (int i = 0; i < nrec; i++) + { + wavebndline.time = (Param.rtlength)*i; + + //slbndline = readBSHline(line); + wavebndvec.push_back(wavebndline); + } + + + return wavebndvec; +} + + +void read_reuse_bndnc(XBGPUParam Param, int rec, float &Trep, double * &qfile, double * &Stfile) +{ + //read the dimentions of grid, levels and time + int status; + int ncid, theta_dimid, tt_dimid, yy_dimid, ee_varid,qx_varid, qy_varid; + + static size_t nrec, ntheta, ntt,nyy; + status = nc_open(Param.wavebndfile.c_str(), NC_NOWRITE, &ncid); + + + status = nc_inq_dimid(ncid, "theta", &theta_dimid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimlen(ncid, theta_dimid, &ntheta); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimid(ncid, "tt", &tt_dimid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimlen(ncid, tt_dimid, &ntt); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimid(ncid, "yy", &yy_dimid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_dimlen(ncid, yy_dimid, &nyy); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_varid(ncid, "ee_bnd", &ee_varid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_varid(ncid, "qx_bnd", &qx_varid); + if (status != NC_NOERR) handle_error(status); + + status = nc_inq_varid(ncid, "qy_bnd", &qy_varid); + if (status != NC_NOERR) handle_error(status); + + int tslenbc = ceil(Param.rtlength / Param.dtbc); + + //Sanity check + if (nyy > Param.ny) + { + // Error + //exit + } + if (ntt > tslenbc) + { + // Error + //exit + } + + double * qxtemp, *qytemp, *eetemp; + //Allocate the temporary array + qxtemp = (double *)malloc(nyy*ntt*sizeof(double)); + qytemp = (double *)malloc(nyy*ntt*sizeof(double)); + eetemp = (double *)malloc(nyy*ntheta*ntt*sizeof(double)); + + static size_t start_ee[] = { rec, 0, 0, 0 }; // start at first value + static size_t count_ee[] = { 1, ntheta, nyy, ntt }; + static size_t start_q[] = { rec, 0, 0 }; // start at first value + static size_t count_q[] = { 1, nyy, ntt }; + + status = nc_get_vara_double(ncid, ee_varid, start_ee, count_ee, eetemp); + if (status != NC_NOERR) handle_error(status); + + status = nc_get_vara_double(ncid, qx_varid, start_q, count_q, qxtemp); + if (status != NC_NOERR) handle_error(status); + + status = nc_get_vara_double(ncid, qy_varid, start_q, count_q, qytemp); + if (status != NC_NOERR) handle_error(status); + + + //close file + status = nc_close(ncid); + + for (int j = 0; j < nyy; j++) + { + for (int m = 0; m < ntt; m++) + { + qfile[j + 0 * nyy + m*nyy * 4]=qxtemp[m + j*ntt] ; + qfile[j + 1 * nyy + m*nyy * 4]=qytemp[m + j*ntt] ; + + for (int itheta = 0; itheta < Param.ntheta; itheta++) + { + Stfile[j + itheta*nyy + m*nyy*Param.ntheta] = eetemp[m + j*ntt + itheta*nyy*ntt]; + } + } + } + + Trep = 15.0; /// Temporary debug + + free(qxtemp); + free(qytemp); + free(eetemp); + } \ No newline at end of file