diff --git a/src/Arrays.h b/src/Arrays.h index 65e9ca35..954ba225 100755 --- a/src/Arrays.h +++ b/src/Arrays.h @@ -59,6 +59,15 @@ struct AdvanceP }; +struct outP +{ + float* z; + short* z_s; + int level; + double xmin, xmax, ymin, ymax; +}; + + struct maskinfo { diff --git a/src/Param.h b/src/Param.h index 734c9f00..a49abdad 100644 --- a/src/Param.h +++ b/src/Param.h @@ -207,6 +207,9 @@ class Param { std::string reftime = ""; // Reference time string as yyyy-mm-ddTHH:MM:SS std::string crs_ref = "no_crs"; //"PROJCS[\"NZGD2000 / New Zealand Transverse Mercator 2000\",GEOGCS[\"NZGD2000\",DATUM[\"New_Zealand_Geodetic_Datum_2000\",SPHEROID[\"GRS 1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4167\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",173],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",1600000],PARAMETER[\"false_northing\",10000000],UNIT[\"metre\",1],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"2193\"]]"; + + bool savebyblk=true; + }; diff --git a/src/Poly.cu b/src/Poly.cu index c189c4f8..53665a79 100644 --- a/src/Poly.cu +++ b/src/Poly.cu @@ -212,7 +212,7 @@ Polygon CounterCWPoly(Polygon Poly) if (sum > 0.0) { log(" Reversing Polygon handedness"); - for (int i = Poly.vertices.size(); i > 0; i--) + for (int i = Poly.vertices.size()-1; i > 0; i--) { // diff --git a/src/ReadInput.cu b/src/ReadInput.cu index 72cd393e..3d9f0bf5 100644 --- a/src/ReadInput.cu +++ b/src/ReadInput.cu @@ -736,6 +736,13 @@ Param readparamstr(std::string line, Param param) param.crs_ref = parametervalue; } + paramvec = { "savebyblk", "writebyblk","saveperblk", "writeperblk","savebyblock", "writebyblock","saveperblock", "writeperblock" }; + parametervalue = findparameter(paramvec, line); + if (!parametervalue.empty()) + { + param.savebyblk = readparambool(parametervalue, param.savebyblk); + } + return param; } diff --git a/src/Updateforcing.cu b/src/Updateforcing.cu index 4a9c9b8c..efb0d497 100755 --- a/src/Updateforcing.cu +++ b/src/Updateforcing.cu @@ -746,11 +746,11 @@ template void deformstep(Param XParam, Loop XLoop, std::vector>> (XParam, XModel.blocks, deform[nd], scale, XModel.evolv.zs, XModel.zb); + AddDeformGPU <<>> (XParam, XModel.blocks, deform[nd], XModel.evolv, scale, XModel.zb); CUDA_CHECK(cudaDeviceSynchronize()); } @@ -779,7 +779,7 @@ template void deformstep(Param XParam, Loop XLoop, std::vector __global__ void AddDeformGPU(Param XParam, BlockP XBlock, deformmap defmap, T scale, T* zs, T* zb) +template __global__ void AddDeformGPU(Param XParam, BlockP XBlock, deformmap defmap, EvolvingP XEv, T scale, T* zb) { unsigned int ix = threadIdx.x; unsigned int iy = threadIdx.y; @@ -801,7 +801,7 @@ template __global__ void AddDeformGPU(Param XParam, BlockP XBlock, // printf("x=%f, y=%f, def=%f\n ", x, y, def); //} - zss = zs[i] + def * scale; + zss = XEv.zs[i] + def * scale; if (defmap.iscavity == true) { zbb = min(zss, zb[i]); @@ -811,7 +811,8 @@ template __global__ void AddDeformGPU(Param XParam, BlockP XBlock, zbb = zb[i] + def * scale; } - zs[i] = zss; + XEv.h[i] = zss - zbb; + XEv.zs[i] = zss; zb[i] = zbb; //zs[i] = zs[i] + def * scale; @@ -821,7 +822,7 @@ template __global__ void AddDeformGPU(Param XParam, BlockP XBlock, } -template __host__ void AddDeformCPU(Param XParam, BlockP XBlock, deformmap defmap, T scale, T* zs, T* zb) +template __host__ void AddDeformCPU(Param XParam, BlockP XBlock, deformmap defmap, EvolvingP XEv, T scale, T* zb) { int ib; @@ -846,7 +847,7 @@ template __host__ void AddDeformCPU(Param XParam, BlockP XBlock, de def = interp2BUQ(x, y, defmap); - zss = zs[i] + def * scale; + zss = XEv.zs[i] + def * scale; if (defmap.iscavity == true) { zbb = min(zss, zb[i]); @@ -856,7 +857,8 @@ template __host__ void AddDeformCPU(Param XParam, BlockP XBlock, de zbb = zb[i] + def * scale; } - zs[i] = zss; + XEv.zs[i] = zss; + XEv.h[i] = zss - zbb; zb[i] = zbb; } } diff --git a/src/Updateforcing.h b/src/Updateforcing.h index 313652d8..a5a8b912 100755 --- a/src/Updateforcing.h +++ b/src/Updateforcing.h @@ -36,7 +36,7 @@ template __host__ void AddRiverForcing(Param XParam, Loop XLoop, st template void deformstep(Param XParam, Loop XLoop, std::vector> deform, Model XModel, Model XModel_g); template __global__ void InjectRiverGPU(Param XParam, River XRiver, T qnow, int* Riverblks, BlockP XBlock, AdvanceP XAdv); -template __global__ void AddDeformGPU(Param XParam, BlockP XBlock, deformmap defmap, T scale, T* zs, T* zb); +template __global__ void AddDeformGPU(Param XParam, BlockP XBlock, deformmap defmap, EvolvingP XEv, T scale, T* zb); #endif diff --git a/src/Write_netcdf.cu b/src/Write_netcdf.cu index 1a1a55e2..b5fe545f 100644 --- a/src/Write_netcdf.cu +++ b/src/Write_netcdf.cu @@ -615,7 +615,7 @@ template void defncvarBUQ(Param XParam, int* activeblk, int* level, T* int shuffle = 1; int deflate = 1; // This switches compression on (1) or off (0). - int deflate_level = 9; // This is the compression level in range 1 (less) - 9 (more). + int deflate_level = 5; // This is the compression level in range 1 (less) - 9 (more). nc_def_var_deflate(ncid, var_id, shuffle, deflate, deflate_level); } @@ -745,6 +745,345 @@ template void defncvarBUQ(Param XParam, int* activeblk, int* level, float template void defncvarBUQ(Param XParam, int* activeblk, int* level, double* blockxo, double* blockyo, std::string varst, int vdim, double* var, outzoneB Xzone); +template void defncvarBUQlev(Param XParam, int* activeblk, int* level, T* blockxo, T* blockyo, std::string varst, std::string longname, std::string stdname, std::string unit, int vdim, T* var, outzoneB Xzone) +{ + + int smallnc = XParam.smallnc; + float scalefactor = XParam.scalefactor; + float addoffset = XParam.addoffset; + //int nx = ceil(XParam.nx / 16.0) * 16.0; + //int ny = ceil(XParam.ny / 16.0) * 16.0; + int status; + int ncid, var_id; + int var_dimid2D[2]; + int var_dimid3D[3]; + //int var_dimid4D[4]; + + short* varblk_s; + float* varblk; + int recid, xid, yid; + int bl, ibl; + //size_t ntheta;// nx and ny are stored in XParam not yet for ntheta + + float fillval = 9.9692e+36f; + short fillval_s = (short)round((9.9692e+36f - addoffset) / scalefactor); + //short Sfillval = 32767; + //short fillval = 32767 + static size_t start2D[] = { 0, 0 }; // start at first value + //static size_t count2D[] = { ny, nx }; + static size_t count2D[] = { (size_t)XParam.blkwidth, (size_t)XParam.blkwidth }; + + static size_t start3D[] = { 0, 0, 0 }; // start at first value + //static size_t count3D[] = { 1, ny, nx }; + static size_t count3D[] = { 1, (size_t)XParam.blkwidth, (size_t)XParam.blkwidth }; + //size_t count3D[3]; + //count3D[0] = 1; + //count3D[1] = XParam.blkwidth; + //count3D[2] = XParam.blkwidth; + + //int minlevzone, maxlevzone; + + std::string outfile = Xzone.outname; + std::vector activeblkzone = Calcactiveblockzone(XParam, activeblk, Xzone); + //Calclevelzone(XParam, minlevzone, maxlevzone, Xzone, level); + + + nc_type VarTYPE; + + if (smallnc > 0) + { + VarTYPE = NC_SHORT; + } + else + { + VarTYPE = NC_FLOAT; + } + + //printf("\n ib=%d count3D=[%d,%d,%d]\n", count3D[0], count3D[1], count3D[2]); + + + status = nc_open(outfile.c_str(), NC_WRITE, &ncid); + if (status != NC_NOERR) handle_ncerror(status); + status = nc_redef(ncid); + if (status != NC_NOERR) handle_ncerror(status); + //Inquire dimensions ids + status = nc_inq_unlimdim(ncid, &recid);//time + if (status != NC_NOERR) handle_ncerror(status); + + varblk = (float*)malloc(XParam.blkwidth * XParam.blkwidth * sizeof(float)); + if (smallnc > 0) + { + + varblk_s = (short*)malloc(XParam.blkwidth * XParam.blkwidth * sizeof(short)); + } + + + std::string xxname, yyname, varname, sign; + + //generate a different variable name for each level and add attribute as necessary + for (int lev = Xzone.minlevel; lev <= Xzone.maxlevel; lev++) + { + + //std::string xxname, yyname, sign; + + lev < 0 ? sign = "N" : sign = "P"; + + + xxname = "xx_" + sign + std::to_string(abs(lev)); + yyname = "yy_" + sign + std::to_string(abs(lev)); + + varname = varst + "_" + sign + std::to_string(abs(lev)); + + + //printf("lev=%d; xxname=%s; yyname=%s;\n", lev, xxname.c_str(), yyname.c_str()); + status = nc_inq_dimid(ncid, xxname.c_str(), &xid); + if (status != NC_NOERR) handle_ncerror(status); + status = nc_inq_dimid(ncid, yyname.c_str(), &yid); + if (status != NC_NOERR) handle_ncerror(status); + + + var_dimid2D[0] = yid; + var_dimid2D[1] = xid; + + var_dimid3D[0] = recid; + var_dimid3D[1] = yid; + var_dimid3D[2] = xid; + + if (vdim == 2) + { + status = nc_def_var(ncid, varname.c_str(), VarTYPE, vdim, var_dimid2D, &var_id); + if (status != NC_NOERR) handle_ncerror(status); + } + else if (vdim == 3) + { + status = nc_def_var(ncid, varname.c_str(), VarTYPE, vdim, var_dimid3D, &var_id); + if (status != NC_NOERR) handle_ncerror(status); + } + + if (smallnc > 0) + { + + status = nc_put_att_short(ncid, var_id, "_FillValue", NC_SHORT, 1, &fillval_s); + if (status != NC_NOERR) handle_ncerror(status); + status = nc_put_att_short(ncid, var_id, "missingvalue", NC_SHORT, 1, &fillval_s); + + if (status != NC_NOERR) handle_ncerror(status); + } + else + { + status = nc_put_att_float(ncid, var_id, "_FillValue", NC_FLOAT, 1, &fillval); + if (status != NC_NOERR) handle_ncerror(status); + status = nc_put_att_float(ncid, var_id, "missingvalue", NC_FLOAT, 1, &fillval); + + if (status != NC_NOERR) handle_ncerror(status); + } + + + if (smallnc > 0) + { + + status = nc_put_att_float(ncid, var_id, "scale_factor", NC_FLOAT, 1, &scalefactor); + if (status != NC_NOERR) handle_ncerror(status); + status = nc_put_att_float(ncid, var_id, "add_offset", NC_FLOAT, 1, &addoffset); + if (status != NC_NOERR) handle_ncerror(status); + } + + + status = nc_put_att_text(ncid, var_id, "standard_name", stdname.size(), stdname.c_str()); + status = nc_put_att_text(ncid, var_id, "long_name", longname.size(), longname.c_str()); + status = nc_put_att_text(ncid, var_id, "units", unit.size(), unit.c_str()); + + std::string crsstrname = "crs"; + status = nc_put_att_text(ncid, var_id, "grid_mapping", crsstrname.size(), crsstrname.c_str()); + + + int shuffle = 1; + int deflate = 1; // This switches compression on (1) or off (0). + int deflate_level = 5; // This is the compression level in range 1 (less) - 9 (more). + nc_def_var_deflate(ncid, var_id, shuffle, deflate, deflate_level); + + } + // End definition + status = nc_enddef(ncid); + if (status != NC_NOERR) handle_ncerror(status); + + //printf("\n ib=%d count3D=[%d,%d,%d]\n", count3D[0], count3D[1], count3D[2]); + + // Now write the initial value of the Variable out + + //std::vector activeblkzone = Calcactiveblockzone(XParam, activeblk, Xzone); + + //#################### + // Create empty array for each levels + std::vector varlayers; + int nx, ny; + for (int levi = Xzone.minlevel; levi <= Xzone.maxlevel; levi++) + { + varlayers.push_back(outP()); + int levindex = (levi - Xzone.minlevel); + + Calcnxnyzone(XParam, levi, nx, ny, Xzone); + varlayers[levindex].z = (float*)malloc(nx * ny * sizeof(float)); + if (smallnc > 0) + { + + varlayers[levindex].z_s = (short*)malloc(nx * ny * sizeof(short)); + } + varlayers[levindex].level = levi; + + for (int j = 0; j < ny; j++) + { + for (int i = 0; i < nx; i++) + { + int n = i + j * nx; + + varlayers[levindex].z[n] = fillval; + if (smallnc > 0) + { + varlayers[levindex].z_s[n] = fillval_s; + } + } + } + } + + //std::string xxname, yyname, varname, sign; + //std::vector activeblkzone = Calcactiveblockzone(XParam, activeblk, Xzone); + + + //int lev, bl; + for (int ibl = 0; ibl < Xzone.nblk; ibl++) + { + //bl = activeblk[Xzone.blk[ibl]]; + //for (int ibl = 0; ibl < XParam.nblk; ibl++) + //{ + bl = activeblkzone[ibl]; + int lev = level[bl]; + + double xxmin, yymin; + int nxlev, nylev; + //double xxmax, yymax; + double initdx = calcres(XParam.dx, XParam.initlevel); + + int io, jo; + //xxmax = Xzone.xmax - calcres(XParam.dx, lev) / 2.0; + //yymax = Xzone.ymax - calcres(XParam.dx, lev) / 2.0; + + int levindex = (lev - Xzone.minlevel); + + Calcnxnyzone(XParam, lev, nxlev, nylev, Xzone); + xxmin = Xzone.xo + calcres(XParam.dx, lev) / 2.0; + yymin = Xzone.yo + calcres(XParam.dx, lev) / 2.0; + + + + jo = round((XParam.yo + blockyo[bl] - yymin) / calcres(XParam.dx, lev)); + io = round((XParam.xo + blockxo[bl] - xxmin) / calcres(XParam.dx, lev)); + + + for (int j = 0; j < XParam.blkwidth; j++) + { + for (int i = 0; i < XParam.blkwidth; i++) + { + int n = (i + XParam.halowidth + XParam.outishift) + (j + XParam.halowidth + XParam.outjshift) * XParam.blkmemwidth + bl * XParam.blksize; + int r = (io + i) + (jo + j) * nxlev; + if (smallnc > 0) + { + // packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) + varlayers[levindex].z_s[r] = (short)round((var[n] - addoffset) / scalefactor); + } + else + { + varlayers[levindex].z[r] = (float)var[n]; + } + } + } + + } + + + for (int levi = Xzone.minlevel; levi <= Xzone.maxlevel; levi++) + { + int nxlev, nylev; + Calcnxnyzone(XParam, levi, nxlev, nylev, Xzone); + //double xxmin, yymin; + levi < 0 ? sign = "N" : sign = "P"; + varname = varst + "_" + sign + std::to_string(abs(levi)); + + status = nc_inq_varid(ncid, varname.c_str(), &var_id); + if (status != NC_NOERR) handle_ncerror(status); + + //xxmin = Xzone.xo + calcres(XParam.dx, levi) / 2.0; + //yymin = Xzone.yo + calcres(XParam.dx, levi) / 2.0; + + int levindex = (levi - Xzone.minlevel); + + if (vdim == 2) + { + + + start2D[0] = 0; // (size_t)round((XParam.yo - yymin) / calcres(XParam.dx, levi)); + start2D[1] = 0; // (size_t)round((XParam.xo - xxmin) / calcres(XParam.dx, levi)); + + count2D[0] = (size_t)nylev; + count2D[1] = (size_t)nxlev; + + if (smallnc > 0) + { + + status = nc_put_vara_short(ncid, var_id, start2D, count2D, varlayers[levindex].z_s); + if (status != NC_NOERR) handle_ncerror(status); + } + else + { + status = nc_put_vara_float(ncid, var_id, start2D, count2D, varlayers[levindex].z); + if (status != NC_NOERR) handle_ncerror(status); + } + } + else if (vdim == 3) + { + start3D[1] = 0;// (size_t)round((XParam.yo - yymin) / calcres(XParam.dx, levi)); + start3D[2] = 0;// (size_t)round((XParam.xo - xxmin) / calcres(XParam.dx, levi)); + + count3D[1] = (size_t)nylev; + count3D[2] = (size_t)nxlev; + + if (smallnc > 0) + { + status = nc_put_vara_short(ncid, var_id, start3D, count3D, varlayers[levindex].z_s); + if (status != NC_NOERR) handle_ncerror(status); + } + else + { + status = nc_put_vara_float(ncid, var_id, start3D, count3D, varlayers[levindex].z); + if (status != NC_NOERR) handle_ncerror(status); + //printf("\n ib=%d start=[%d,%d,%d]; initlevel=%d; initdx=%f; level=%d; xo=%f; yo=%f; blockxo[ib]=%f xxmin=%f blockyo[ib]=%f yymin=%f startfl=%f\n", bl, start3D[0], start3D[1], start3D[2], XParam.initlevel, initdx, lev, Xzone.xo, Xzone.yo, blockxo[bl], xxmin, blockyo[bl], yymin, (blockyo[bl] - yymin) / calcres(XParam.dx, lev)); + //printf("\n varblk[0]=%f varblk[255]=%f\n", varblk[0], varblk[255]); + //printf("\n ib=%d count3D=[%d,%d,%d]\n", count3D[0], count3D[1], count3D[2]); + //printf("\n ib=%d; level=%d; blockxo[ib]=%f blockyo[ib]=%f \n", bl, lev, blockxo[bl], blockyo[bl]); + } + + } + + } + + for (int levi = Xzone.minlevel; levi <= Xzone.maxlevel; levi++) + { + int levindex = (levi - Xzone.minlevel); + if (smallnc > 0) + { + + free(varlayers[levindex].z_s); + } + free(varlayers[levindex].z); + } + //close and save new file + status = nc_close(ncid); + if (status != NC_NOERR) handle_ncerror(status); + +} + + + template void writencvarstepBUQ(Param XParam, int vdim, int * activeblk, int* level, T * blockxo, T *blockyo, std::string varst, T * var, outzoneB Xzone) { @@ -915,6 +1254,219 @@ extern "C" void writenctimestep(std::string outfile, double totaltime) if (status != NC_NOERR) handle_ncerror(status); } + +template void writencvarstepBUQlev(Param XParam, int vdim, int* activeblk, int* level, T* blockxo, T* blockyo, std::string varst, T* var, outzoneB Xzone) +{ + int status, ncid, recid, var_id; + static size_t nrec; + short* varblk_s; + float* varblk; + //int nx, ny; + //int dimids[NC_MAX_VAR_DIMS]; + //size_t *ddim, *start, *count; + //XParam.outfile.c_str() + + + float fillval = 9.9692e+36f; + short fillval_s = (short)round((9.9692e+36f - XParam.addoffset) / XParam.scalefactor); + + static size_t start2D[] = { 0, 0 }; // start at first value + //static size_t count2D[] = { ny, nx }; + static size_t count2D[] = { (size_t)XParam.blkwidth, (size_t)XParam.blkwidth }; + + static size_t start3D[] = { 0, 0, 0 }; // start at first value // This is updated to nrec-1 further down + //static size_t count3D[] = { 1, ny, nx }; + static size_t count3D[] = { 1, (size_t)XParam.blkwidth, (size_t)XParam.blkwidth }; + + int smallnc = XParam.smallnc; + float scalefactor = XParam.scalefactor; + float addoffset = XParam.addoffset; + + status = nc_open(Xzone.outname.c_str(), NC_WRITE, &ncid); + if (status != NC_NOERR) handle_ncerror(status); + //read id from time dimension + status = nc_inq_unlimdim(ncid, &recid); + if (status != NC_NOERR) handle_ncerror(status); + status = nc_inq_dimlen(ncid, recid, &nrec); + if (status != NC_NOERR) handle_ncerror(status); + + start3D[0] = nrec - 1; + + + // Create empty array for each levels + std::vector varlayers; + int nx, ny; + for (int levi = Xzone.minlevel; levi <= Xzone.maxlevel; levi++) + { + varlayers.push_back(outP()); + int levindex = (levi - Xzone.minlevel); + + Calcnxnyzone(XParam, levi, nx, ny, Xzone); + varlayers[levindex].z = (float*)malloc(nx * ny * sizeof(float)); + if (smallnc > 0) + { + + varlayers[levindex].z_s = (short*)malloc(nx * ny * sizeof(short)); + } + varlayers[levindex].level = levi; + + for (int j = 0; j < ny; j++) + { + for (int i = 0; i < nx; i++) + { + int n = i + j * nx; + + varlayers[levindex].z[n] = fillval; + if (smallnc > 0) + { + varlayers[levindex].z_s[n] = fillval_s; + } + } + } + } + + std::string xxname, yyname, varname, sign; + std::vector activeblkzone = Calcactiveblockzone(XParam, activeblk, Xzone); + + + int lev, bl; + for (int ibl = 0; ibl < Xzone.nblk; ibl++) + { + //bl = activeblk[Xzone.blk[ibl]]; + //for (int ibl = 0; ibl < XParam.nblk; ibl++) + //{ + bl = activeblkzone[ibl]; + lev = level[bl]; + + double xxmin, yymin; + int nxlev, nylev; + //double xxmax, yymax; + double initdx = calcres(XParam.dx, XParam.initlevel); + + int io, jo; + //xxmax = Xzone.xmax - calcres(XParam.dx, lev) / 2.0; + //yymax = Xzone.ymax - calcres(XParam.dx, lev) / 2.0; + + int levindex = (lev - Xzone.minlevel); + + Calcnxnyzone(XParam, lev, nxlev, nylev, Xzone); + xxmin = Xzone.xo + calcres(XParam.dx, lev) / 2.0; + yymin = Xzone.yo + calcres(XParam.dx, lev) / 2.0; + + + + jo = round((XParam.yo + blockyo[bl] - yymin) / calcres(XParam.dx, lev)); + io = round((XParam.xo + blockxo[bl] - xxmin) / calcres(XParam.dx, lev)); + + + for (int j = 0; j < XParam.blkwidth; j++) + { + for (int i = 0; i < XParam.blkwidth; i++) + { + int n = (i + XParam.halowidth + XParam.outishift) + (j + XParam.halowidth + XParam.outjshift) * XParam.blkmemwidth + bl * XParam.blksize; + int r = (io + i) + (jo + j) * nxlev; + if (smallnc > 0) + { + // packed_data_value = nint((unpacked_data_value - add_offset) / scale_factor) + varlayers[levindex].z_s[r] = (short)round((var[n] - addoffset) / scalefactor); + } + else + { + varlayers[levindex].z[r] = (float)var[n]; + } + } + } + + } + + + for (int levi = Xzone.minlevel; levi <= Xzone.maxlevel; levi++) + { + int nxlev, nylev; + Calcnxnyzone(XParam, levi, nxlev, nylev, Xzone); + //double xxmin, yymin; + levi < 0 ? sign = "N" : sign = "P"; + varname = varst + "_" + sign + std::to_string(abs(levi)); + + status = nc_inq_varid(ncid, varname.c_str(), &var_id); + if (status != NC_NOERR) handle_ncerror(status); + + //xxmin = Xzone.xo + calcres(XParam.dx, levi) / 2.0; + //yymin = Xzone.yo + calcres(XParam.dx, levi) / 2.0; + + int levindex = (levi - Xzone.minlevel); + + if (vdim == 2) + { + + + start2D[0] = 0; // (size_t)round((XParam.yo - yymin) / calcres(XParam.dx, levi)); + start2D[1] = 0; // (size_t)round((XParam.xo - xxmin) / calcres(XParam.dx, levi)); + + count2D[0] = (size_t)nylev; + count2D[1] = (size_t)nxlev; + + if (smallnc > 0) + { + + status = nc_put_vara_short(ncid, var_id, start2D, count2D, varlayers[levindex].z_s); + if (status != NC_NOERR) handle_ncerror(status); + } + else + { + status = nc_put_vara_float(ncid, var_id, start2D, count2D, varlayers[levindex].z); + if (status != NC_NOERR) handle_ncerror(status); + } + } + else if (vdim == 3) + { + start3D[1] = 0;// (size_t)round((XParam.yo - yymin) / calcres(XParam.dx, levi)); + start3D[2] = 0;// (size_t)round((XParam.xo - xxmin) / calcres(XParam.dx, levi)); + + count3D[1] = (size_t)nylev; + count3D[2] = (size_t)nxlev; + + if (smallnc > 0) + { + status = nc_put_vara_short(ncid, var_id, start3D, count3D, varlayers[levindex].z_s); + if (status != NC_NOERR) handle_ncerror(status); + } + else + { + status = nc_put_vara_float(ncid, var_id, start3D, count3D, varlayers[levindex].z); + if (status != NC_NOERR) handle_ncerror(status); + //printf("\n ib=%d start=[%d,%d,%d]; initlevel=%d; initdx=%f; level=%d; xo=%f; yo=%f; blockxo[ib]=%f xxmin=%f blockyo[ib]=%f yymin=%f startfl=%f\n", bl, start3D[0], start3D[1], start3D[2], XParam.initlevel, initdx, lev, Xzone.xo, Xzone.yo, blockxo[bl], xxmin, blockyo[bl], yymin, (blockyo[bl] - yymin) / calcres(XParam.dx, lev)); + //printf("\n varblk[0]=%f varblk[255]=%f\n", varblk[0], varblk[255]); + //printf("\n ib=%d count3D=[%d,%d,%d]\n", count3D[0], count3D[1], count3D[2]); + //printf("\n ib=%d; level=%d; blockxo[ib]=%f blockyo[ib]=%f \n", bl, lev, blockxo[bl], blockyo[bl]); + } + + } + + } + + for (int levi = Xzone.minlevel; levi <= Xzone.maxlevel; levi++) + { + int levindex = (levi - Xzone.minlevel); + if (smallnc > 0) + { + + free(varlayers[levindex].z_s); + } + free(varlayers[levindex].z); + } + //close and save new file + status = nc_close(ncid); + if (status != NC_NOERR) handle_ncerror(status); +} + +// Scope for compiler to know what function to compile + +template void writencvarstepBUQlev(Param XParam, int vdim, int* activeblk, int* level, float* blockxo, float* blockyo, std::string varst, float* var, outzoneB Xzone); +template void writencvarstepBUQlev(Param XParam, int vdim, int* activeblk, int* level, double* blockxo, double* blockyo, std::string varst, double* var, outzoneB Xzone); + + + template void InitSave2Netcdf(Param &XParam, Model &XModel) { if (!XParam.outvars.empty()) @@ -928,7 +1480,16 @@ template void InitSave2Netcdf(Param &XParam, Model &XModel) for (int ivar = 0; ivar < XParam.outvars.size(); ivar++) { std::string varstr = XParam.outvars[ivar]; - defncvarBUQ(XParam, XModel.blocks.active, XModel.blocks.level, XModel.blocks.xo, XModel.blocks.yo, varstr,XModel.Outvarlongname[varstr],XModel.Outvarstdname[varstr],XModel.Outvarunits[varstr], 3, XModel.OutputVarMap[varstr], XModel.blocks.outZone[o]); + if (XParam.savebyblk) + { + defncvarBUQ(XParam, XModel.blocks.active, XModel.blocks.level, XModel.blocks.xo, XModel.blocks.yo, varstr, XModel.Outvarlongname[varstr], XModel.Outvarstdname[varstr], XModel.Outvarunits[varstr], 3, XModel.OutputVarMap[varstr], XModel.blocks.outZone[o]); + + } + else + { + defncvarBUQlev(XParam, XModel.blocks.active, XModel.blocks.level, XModel.blocks.xo, XModel.blocks.yo, varstr, XModel.Outvarlongname[varstr], XModel.Outvarstdname[varstr], XModel.Outvarunits[varstr], 3, XModel.OutputVarMap[varstr], XModel.blocks.outZone[o]); + + } } } } @@ -947,7 +1508,14 @@ template void Save2Netcdf(Param XParam,Loop XLoop, Model XModel) writenctimestep(XModel.blocks.outZone[o].outname, XLoop.totaltime); for (int ivar = 0; ivar < XParam.outvars.size(); ivar++) { - writencvarstepBUQ(XParam, 3, XModel.blocks.active, XModel.blocks.level, XModel.blocks.xo, XModel.blocks.yo, XParam.outvars[ivar], XModel.OutputVarMap[XParam.outvars[ivar]], XModel.blocks.outZone[o]); + if (XParam.savebyblk) + { + writencvarstepBUQ(XParam, 3, XModel.blocks.active, XModel.blocks.level, XModel.blocks.xo, XModel.blocks.yo, XParam.outvars[ivar], XModel.OutputVarMap[XParam.outvars[ivar]], XModel.blocks.outZone[o]); + } + else + { + writencvarstepBUQlev(XParam, 3, XModel.blocks.active, XModel.blocks.level, XModel.blocks.xo, XModel.blocks.yo, XParam.outvars[ivar], XModel.OutputVarMap[XParam.outvars[ivar]], XModel.blocks.outZone[o]); + } } } } diff --git a/src/Write_netcdf.h b/src/Write_netcdf.h index 717e12b3..0f6e1927 100755 --- a/src/Write_netcdf.h +++ b/src/Write_netcdf.h @@ -8,6 +8,7 @@ #include "ReadInput.h" #include "MemManagement.h" #include "Util_CPU.h" +#include "Arrays.h" void handle_ncerror(int status); template void creatncfileBUQ(Param &XParam, int* activeblk, int* level, T* blockxo, T* blockyo, outzoneB &Xzone);