From 694b1ba32addc732bd58df890a69c57b27e0dd4a Mon Sep 17 00:00:00 2001 From: cenamiller Date: Thu, 15 Aug 2024 11:39:49 -0600 Subject: [PATCH 1/3] Added outputAnalysis_thermo to CostFunctionXYZ.cpp and added in mode signifier (wind or thermo) to CostFunctionXYZ::outputAnalysis_thermo, CostFunctionXYZ::outputAnalysis, and CostFuctionRTZ::outputAnalysis --- src/CostFunctionRTZ.cpp | 7 +- src/CostFunctionXYZ.cpp | 247 +++++++++++++++++++++++++++++++++++++++- src/CostFunctionXYZ.h | 1 + 3 files changed, 249 insertions(+), 6 deletions(-) diff --git a/src/CostFunctionRTZ.cpp b/src/CostFunctionRTZ.cpp index 0914d94..ae3a396 100644 --- a/src/CostFunctionRTZ.cpp +++ b/src/CostFunctionRTZ.cpp @@ -24,9 +24,10 @@ CostFunctionRTZ::~CostFunctionRTZ() } bool CostFunctionRTZ::outputAnalysis(const std::string& suffix, real* Astate) { + std::string samuraiMode = "wind"; cout << "Outputting " << suffix << "...\n"; // H --> to Mish for output - std::string samuraiout = "samurai_RTZ_" + suffix + ".out"; + std::string samuraiout = "samurai_RTZ_" + samuraiMode + "_" + suffix + ".out"; ofstream samuraistream; if ((*configHash)["output_txt"] == "true") { samuraistream.open(outputPath + "/" + samuraiout); @@ -447,12 +448,12 @@ bool CostFunctionRTZ::outputAnalysis(const std::string& suffix, real* Astate) } } - std::string fileName = "samurai_RTZ_" + suffix; + std::string fileName = "samurai_RTZ_" + samuraiMode + "_" + suffix; std::string outFileName = outputPath + "/" + fileName; // Write the Obs to a summary text file if ((*configHash)["output_qc"] == "true") { - std::string qcout = "samurai_QC_" + suffix + ".out"; + std::string qcout = "samurai_QC_" + samuraiMode + "_" + suffix + ".out"; std::string qcFileName = outputPath + "/" + qcout; ofstream qcstream(qcFileName); ostream_iterator os(qcstream, "\t "); diff --git a/src/CostFunctionXYZ.cpp b/src/CostFunctionXYZ.cpp index 6792e8a..26c61c5 100644 --- a/src/CostFunctionXYZ.cpp +++ b/src/CostFunctionXYZ.cpp @@ -507,6 +507,7 @@ bool CostFunctionXYZ::outputAnalysis(const std::string& suffix, real* Astate) bool debug_bgState = isTrue("debug_bgState"); fractl_mode = isEqual("bkgd_obs_interpolation", "fractl"); + std::string samuraiMode = "wind"; if ( debug_bgState) { std::cout << "---- start of debug_bgState" << std::endl; @@ -520,7 +521,7 @@ bool CostFunctionXYZ::outputAnalysis(const std::string& suffix, real* Astate) cout << "Outputting " << suffix << "...\n"; // H --> to Mish for output - std::string samuraiout = "samurai_XYZ_" + suffix + ".out"; + std::string samuraiout = "samurai_XYZ_" + samuraiMode + "_" + suffix + ".out"; ofstream samuraistream; ofstream *samStreamPtr = NULL;; @@ -573,12 +574,12 @@ bool CostFunctionXYZ::outputAnalysis(const std::string& suffix, real* Astate) std::cout << "--- end of debug_ref_state" << std::endl; } - std::string fileName = "samurai_XYZ_" + suffix; + std::string fileName = "samurai_XYZ_" + samuraiMode + "_" + suffix; std::string outFileName = outputPath + "/" + fileName; // Write the Obs to a summary text file if ((*configHash)["output_qc"] == "true") { - std::string qcout = "samurai_QC_" + suffix + ".out"; + std::string qcout = "samurai_QC_" + samuraiMode + "_" + suffix + ".out"; std::string qcFileName = outputPath + "/" + qcout; ofstream qcstream(qcFileName); ostream_iterator os(qcstream, "\t "); @@ -695,6 +696,246 @@ bool CostFunctionXYZ::outputAnalysis(const std::string& suffix, real* Astate) } +bool CostFunctionXYZ::outputAnalysis_thermo(const std::string& suffix, real* Astate) +{ + std::string samuraiMode = "thermo"; + cout << "Outputting " << suffix.toStdString() << "...\n"; + // H --> to Mish for output + QString samuraiout = "samurai_XYZ_" + samuraiMode + "_" + suffix + ".out"; + ofstream samuraistream; + if (configHash->value("output_txt") == "true") { + samuraistream.open(outputPath.absoluteFilePath(samuraiout).toAscii().data()); + samuraistream << "X\tY\tZ\tpip\tthetarhop\tftheta\n"; + samuraistream.precision(10); + } + + int analysisDim = 12; + int analysisSize = (iDim-2)*(jDim-2)*(kDim-2); + finalAnalysis = new real[analysisSize*analysisDim]; + real gausspoint = 0.5*sqrt(1./3.); + + for (int iIndex = 1; iIndex < iDim-1; iIndex++) { + for (int ihalf = 0; ihalf <= mishFlag; ihalf++) { + for (int imu = -ihalf; imu <= ihalf; imu++) { + real i = iMin + DI * (iIndex + (gausspoint * imu + 0.5*ihalf)); + if (i > ((iDim-1)*DI + iMin)) continue; + + for (int jIndex = 1; jIndex < jDim-1; jIndex++) { + for (int jhalf =0; jhalf <=mishFlag; jhalf++) { + for (int jmu = -jhalf; jmu <= jhalf; jmu++) { + real j = jMin + DJ * (jIndex + (gausspoint * jmu + 0.5*jhalf)); + if (j > ((jDim-1)*DJ + jMin)) continue; + + for (int kIndex = 1; kIndex < kDim-1; kIndex++) { + for (int khalf =0; khalf <=mishFlag; khalf++) { + for (int kmu = -khalf; kmu <= khalf; kmu++) { + real k = kMin + DK * (kIndex + (gausspoint * kmu + 0.5*khalf)); + if (k > ((kDim-1)*DK + kMin)) continue; + + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0.; + real jbasis = 0.; + real kbasis = 0.; + real idbasis = 0.; + real jdbasis = 0.; + real kdbasis = 0.; + real pip = 0.; + real thetarhop = 0.; + real ftheta = 0.; + real pipdx = 0.; real thetarhopdx = 0.; real fthetadx = 0.; + real pipdy = 0.; real thetarhopdy = 0.; real fthetady = 0.; + real pipdz = 0.; real thetarhopdz = 0.; real fthetadz = 0.; + + for (int var = 0; var < varDim; var++) { + for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 0, iBCL[var], iBCR[var]); + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 0, jBCL[var], jBCR[var]); + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 0, kBCL[var], kBCR[var]); + idbasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 1, iBCL[var], iBCR[var]); + jdbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 1, jBCL[var], jBCR[var]); + kdbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 1, kBCL[var], kBCR[var]); + real basis3x = ibasis*jbasis*kbasis; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + switch (var) { + case 0: + pip += Astate[aIndex] * basis3x; + pipdx += Astate[aIndex] * idbasis * jbasis * kbasis; + pipdy += Astate[aIndex] * ibasis * jdbasis * kbasis; + pipdz += Astate[aIndex] * ibasis * jbasis * kdbasis; + break; + case 1: + thetarhop += Astate[aIndex + 1] * basis3x; + thetarhopdx += Astate[aIndex + 1] * idbasis * jbasis * kbasis; + thetarhopdy += Astate[aIndex + 1] * ibasis * jdbasis * kbasis; + thetarhopdz += Astate[aIndex + 1] * ibasis * jbasis * kdbasis; + break; + case 2: + ftheta += Astate[aIndex + 2] * basis3x; + fthetadx += Astate[aIndex + 2] * idbasis * jbasis * kbasis; + fthetady += Astate[aIndex + 2] * ibasis * jdbasis * kbasis; + fthetadz += Astate[aIndex + 2] * ibasis * jbasis * kdbasis; + break; + } + } + } + } + } + + + if (configHash->value("output_txt") == "true") { + samuraistream << scientific << i << "\t" << j << "\t" << k + << "\t" << pip << "\t" << thetarhop << "\t" << ftheta << "\t" << "\n"; + } + + // On the nodes + if (!ihalf and !jhalf and !khalf){ + int fIndex = (iDim-2)*(jDim-2)*(kDim-2); + int posIndex = (iDim-2)*(jDim-2)*(kIndex-1) + (iDim-2)*(jIndex-1) + (iIndex-1); + finalAnalysis[fIndex * 0 + posIndex] = pip; + finalAnalysis[fIndex * 1 + posIndex] = thetarhop; + finalAnalysis[fIndex * 2 + posIndex] = ftheta; + finalAnalysis[fIndex * 3 + posIndex] = pipdx; + finalAnalysis[fIndex * 4 + posIndex] = thetarhopdx; + finalAnalysis[fIndex * 5 + posIndex] = fthetadx; + finalAnalysis[fIndex * 6 + posIndex] = pipdy; + finalAnalysis[fIndex * 7 + posIndex] = thetarhopdy; + finalAnalysis[fIndex * 8 + posIndex] = fthetady; + finalAnalysis[fIndex * 9 + posIndex] = pipdz; + finalAnalysis[fIndex * 10 + posIndex] = thetarhopdz; + finalAnalysis[fIndex * 11 + posIndex] = fthetadz; + } + } + } + } + } + } + } + } + } + } + + QString fileName = "samurai_XYZ_" + samuraiMode + "_" + suffix; + QString outFileName = outputPath.absoluteFilePath(fileName); + + // Write the Obs to a summary text file + if (configHash->value("output_qc") == "true") { + QString qcout = "samurai_QC_" + samuraiMode + "_" + suffix + ".out"; + QString qcFileName = outputPath.absoluteFilePath(qcout); + ofstream qcstream(qcFileName.toAscii().data()); + ostream_iterator os(qcstream, "\t "); + *os++ = "Observation"; + *os++ = "Inverse Error"; + *os++ = "X"; + *os++ = "Y"; + *os++ = "Z"; + *os++ = "Type"; + *os++ = "Time"; + *os++ = "pip"; + *os++ = "thetarhop"; + *os++ = "ftheta"; + *os++ = "Analysis"; + *os++ = "Background"; + qcstream << endl; + qcstream.precision(10); + + ostream_iterator od(qcstream, "\t "); + for (int m = 0; m < mObs; m++) { + int mi = m*(obMetaSize+varDim*derivDim); + real i = obsVector[mi+2]; + real j = obsVector[mi+3]; + real k = obsVector[mi+4]; + real tempsum = 0; + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0; + real jbasis = 0; + real kbasis = 0; + for (int var = 0; var < varDim; var++) { + for (int d = 0; d < derivDim; d++) { + int wgt_index = mi + obMetaSize +d*varDim + var; + if (!obsVector[wgt_index]) continue; + for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, derivative[d][2], kBCL[var], kBCR[var]); + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, derivative[d][1], iBCL[var], iBCR[var]); + + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, derivative[d][0], jBCL[var], jBCR[var]); + tempsum += Astate[aIndex + var] * ibasis * jbasis * kbasis * obsVector[wgt_index]; + } + } + } + } + } + for (int t=0; t < obMetaSize-1; t++) { + *od++ = obsVector[mi+t]; + } + int unixtime = (int)obsVector[mi+6]; + QDateTime obtime; + obtime.setTime_t(unixtime); + obtime.setTimeSpec(Qt::UTC); + QString timestring = obtime.toString("hh:mm:ss.zzz"); + qcstream << timestring.toStdString() << "\t"; + + // Multiply the weight by the ob -- Observations.in has individual weights already + // Only non-derivative for now + for (int t=obMetaSize; tvalue("output_netcdf") == "true") { + QString cdfFileName = outFileName + ".nc"; + if (!writeNetCDF(outputPath.absoluteFilePath(cdfFileName))) + cout << "Error writing netcdf file " << cdfFileName.toStdString() << endl; + } + // Write out to an asi file + if (configHash->value("output_asi") == "true") { + QString asiFileName = outFileName + ".asi"; + if (!writeAsi(outputPath.absoluteFilePath(asiFileName))) + cout << "Error writing asi file " << asiFileName.toStdString() << endl; + } + // Set the domain back + adjustInternalDomain(1); + + // Free the memory for the analysis variables + delete[] finalAnalysis; + + return true; + +} + bool CostFunctionXYZ::writeNetCDF(const std::string& netcdfFileName) { Nc3Error err(Nc3Error::verbose_nonfatal); diff --git a/src/CostFunctionXYZ.h b/src/CostFunctionXYZ.h index 69f1b64..894d72c 100644 --- a/src/CostFunctionXYZ.h +++ b/src/CostFunctionXYZ.h @@ -27,6 +27,7 @@ class CostFunctionXYZ: public CostFunction3D private: bool outputAnalysis(const std::string& suffix, real* Astate); + bool outputAnalysis_thermo(const std::string& suffix, real* Astate); bool writeAsi(const std::string& asiFileName); bool writeNetCDF(const std::string& netcdfFileName); bool SItransform(size_t numVars, double *finalAnalysis, double *mishData, real *Astate, ofstream *outStream); From 2f8b996c78a4a89386f8597f003578f2d143f859 Mon Sep 17 00:00:00 2001 From: cenamiller Date: Mon, 19 Aug 2024 15:52:11 -0600 Subject: [PATCH 2/3] Removed stray QTstring library calls from outputAnalysis_thermo --- src/CostFunctionXYZ.cpp | 435 ++++++++++++++++++++-------------------- 1 file changed, 217 insertions(+), 218 deletions(-) diff --git a/src/CostFunctionXYZ.cpp b/src/CostFunctionXYZ.cpp index 26c61c5..ee632d7 100644 --- a/src/CostFunctionXYZ.cpp +++ b/src/CostFunctionXYZ.cpp @@ -698,242 +698,241 @@ bool CostFunctionXYZ::outputAnalysis(const std::string& suffix, real* Astate) bool CostFunctionXYZ::outputAnalysis_thermo(const std::string& suffix, real* Astate) { - std::string samuraiMode = "thermo"; - cout << "Outputting " << suffix.toStdString() << "...\n"; - // H --> to Mish for output - QString samuraiout = "samurai_XYZ_" + samuraiMode + "_" + suffix + ".out"; - ofstream samuraistream; - if (configHash->value("output_txt") == "true") { - samuraistream.open(outputPath.absoluteFilePath(samuraiout).toAscii().data()); - samuraistream << "X\tY\tZ\tpip\tthetarhop\tftheta\n"; - samuraistream.precision(10); - } + std::string samuraiMode = "thermo"; + std::cout << "Outputting " << suffix << "...\n"; + // H --> to Mish for output + std::string samuraiout = "samurai_XYZ_" + samuraiMode + "_" + suffix + ".out"; + std::ofstream samuraistream; + if ((*configHash)["output_txt"] == "true") { + samuraistream.open(outputPath + samuraiout); + samuraistream << "X\tY\tZ\tpip\tthetarhop\tftheta\n"; + samuraistream.precision(10); + } + + int analysisDim = 12; + int analysisSize = (iDim-2)*(jDim-2)*(kDim-2); + finalAnalysis = new real[analysisSize*analysisDim]; + real gausspoint = 0.5*sqrt(1./3.); - int analysisDim = 12; - int analysisSize = (iDim-2)*(jDim-2)*(kDim-2); - finalAnalysis = new real[analysisSize*analysisDim]; - real gausspoint = 0.5*sqrt(1./3.); - - for (int iIndex = 1; iIndex < iDim-1; iIndex++) { - for (int ihalf = 0; ihalf <= mishFlag; ihalf++) { - for (int imu = -ihalf; imu <= ihalf; imu++) { - real i = iMin + DI * (iIndex + (gausspoint * imu + 0.5*ihalf)); - if (i > ((iDim-1)*DI + iMin)) continue; - - for (int jIndex = 1; jIndex < jDim-1; jIndex++) { - for (int jhalf =0; jhalf <=mishFlag; jhalf++) { - for (int jmu = -jhalf; jmu <= jhalf; jmu++) { - real j = jMin + DJ * (jIndex + (gausspoint * jmu + 0.5*jhalf)); - if (j > ((jDim-1)*DJ + jMin)) continue; - - for (int kIndex = 1; kIndex < kDim-1; kIndex++) { - for (int khalf =0; khalf <=mishFlag; khalf++) { - for (int kmu = -khalf; kmu <= khalf; kmu++) { - real k = kMin + DK * (kIndex + (gausspoint * kmu + 0.5*khalf)); - if (k > ((kDim-1)*DK + kMin)) continue; - - int ii = (int)((i - iMin)*DIrecip); - int jj = (int)((j - jMin)*DJrecip); - int kk = (int)((k - kMin)*DKrecip); - real ibasis = 0.; - real jbasis = 0.; - real kbasis = 0.; - real idbasis = 0.; - real jdbasis = 0.; - real kdbasis = 0.; - real pip = 0.; - real thetarhop = 0.; - real ftheta = 0.; - real pipdx = 0.; real thetarhopdx = 0.; real fthetadx = 0.; - real pipdy = 0.; real thetarhopdy = 0.; real fthetady = 0.; - real pipdz = 0.; real thetarhopdz = 0.; real fthetadz = 0.; - - for (int var = 0; var < varDim; var++) { - for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { - for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { - int iNode = iiNode; - if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; - if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); - if ((iNode < 0) or (iNode >= iDim)) continue; - for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { - int jNode = jjNode; - if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; - if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); - if ((jNode < 0) or (jNode >= jDim)) continue; - ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 0, iBCL[var], iBCR[var]); - jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 0, jBCL[var], jBCR[var]); - kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 0, kBCL[var], kBCR[var]); - idbasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 1, iBCL[var], iBCR[var]); - jdbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 1, jBCL[var], jBCR[var]); - kdbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 1, kBCL[var], kBCR[var]); - real basis3x = ibasis*jbasis*kbasis; - int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; - switch (var) { - case 0: - pip += Astate[aIndex] * basis3x; - pipdx += Astate[aIndex] * idbasis * jbasis * kbasis; - pipdy += Astate[aIndex] * ibasis * jdbasis * kbasis; - pipdz += Astate[aIndex] * ibasis * jbasis * kdbasis; - break; - case 1: - thetarhop += Astate[aIndex + 1] * basis3x; - thetarhopdx += Astate[aIndex + 1] * idbasis * jbasis * kbasis; - thetarhopdy += Astate[aIndex + 1] * ibasis * jdbasis * kbasis; - thetarhopdz += Astate[aIndex + 1] * ibasis * jbasis * kdbasis; - break; - case 2: - ftheta += Astate[aIndex + 2] * basis3x; - fthetadx += Astate[aIndex + 2] * idbasis * jbasis * kbasis; - fthetady += Astate[aIndex + 2] * ibasis * jdbasis * kbasis; - fthetadz += Astate[aIndex + 2] * ibasis * jbasis * kdbasis; - break; - } - } - } - } - } - - - if (configHash->value("output_txt") == "true") { - samuraistream << scientific << i << "\t" << j << "\t" << k - << "\t" << pip << "\t" << thetarhop << "\t" << ftheta << "\t" << "\n"; - } - - // On the nodes - if (!ihalf and !jhalf and !khalf){ - int fIndex = (iDim-2)*(jDim-2)*(kDim-2); - int posIndex = (iDim-2)*(jDim-2)*(kIndex-1) + (iDim-2)*(jIndex-1) + (iIndex-1); - finalAnalysis[fIndex * 0 + posIndex] = pip; - finalAnalysis[fIndex * 1 + posIndex] = thetarhop; - finalAnalysis[fIndex * 2 + posIndex] = ftheta; - finalAnalysis[fIndex * 3 + posIndex] = pipdx; - finalAnalysis[fIndex * 4 + posIndex] = thetarhopdx; - finalAnalysis[fIndex * 5 + posIndex] = fthetadx; - finalAnalysis[fIndex * 6 + posIndex] = pipdy; - finalAnalysis[fIndex * 7 + posIndex] = thetarhopdy; - finalAnalysis[fIndex * 8 + posIndex] = fthetady; - finalAnalysis[fIndex * 9 + posIndex] = pipdz; - finalAnalysis[fIndex * 10 + posIndex] = thetarhopdz; - finalAnalysis[fIndex * 11 + posIndex] = fthetadz; - } - } - } + for (int iIndex = 1; iIndex < iDim-1; iIndex++) { + for (int ihalf = 0; ihalf <= mishFlag; ihalf++) { + for (int imu = -ihalf; imu <= ihalf; imu++) { + real i = iMin + DI * (iIndex + (gausspoint * imu + 0.5*ihalf)); + if (i > ((iDim-1)*DI + iMin)) continue; + + for (int jIndex = 1; jIndex < jDim-1; jIndex++) { + for (int jhalf =0; jhalf <=mishFlag; jhalf++) { + for (int jmu = -jhalf; jmu <= jhalf; jmu++) { + real j = jMin + DJ * (jIndex + (gausspoint * jmu + 0.5*jhalf)); + if (j > ((jDim-1)*DJ + jMin)) continue; + + for (int kIndex = 1; kIndex < kDim-1; kIndex++) { + for (int khalf =0; khalf <=mishFlag; khalf++) { + for (int kmu = -khalf; kmu <= khalf; kmu++) { + real k = kMin + DK * (kIndex + (gausspoint * kmu + 0.5*khalf)); + if (k > ((kDim-1)*DK + kMin)) continue; + + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0.; + real jbasis = 0.; + real kbasis = 0.; + real idbasis = 0.; + real jdbasis = 0.; + real kdbasis = 0.; + real pip = 0.; + real thetarhop = 0.; + real ftheta = 0.; + real pipdx = 0.; real thetarhopdx = 0.; real fthetadx = 0.; + real pipdy = 0.; real thetarhopdy = 0.; real fthetady = 0.; + real pipdz = 0.; real thetarhopdz = 0.; real fthetadz = 0.; + + for (int var = 0; var < varDim; var++) { + for (int kNode = std::max(kk-1,0); kNode <= std::min(kk+2,kDim-1); ++kNode) { + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 0, iBCL[var], iBCR[var]); + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 0, jBCL[var], jBCR[var]); + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 0, kBCL[var], kBCR[var]); + idbasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 1, iBCL[var], iBCR[var]); + jdbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 1, jBCL[var], jBCR[var]); + kdbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 1, kBCL[var], kBCR[var]); + real basis3x = ibasis*jbasis*kbasis; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + switch (var) { + case 0: + pip += Astate[aIndex] * basis3x; + pipdx += Astate[aIndex] * idbasis * jbasis * kbasis; + pipdy += Astate[aIndex] * ibasis * jdbasis * kbasis; + pipdz += Astate[aIndex] * ibasis * jbasis * kdbasis; + break; + case 1: + thetarhop += Astate[aIndex + 1] * basis3x; + thetarhopdx += Astate[aIndex + 1] * idbasis * jbasis * kbasis; + thetarhopdy += Astate[aIndex + 1] * ibasis * jdbasis * kbasis; + thetarhopdz += Astate[aIndex + 1] * ibasis * jbasis * kdbasis; + break; + case 2: + ftheta += Astate[aIndex + 2] * basis3x; + fthetadx += Astate[aIndex + 2] * idbasis * jbasis * kbasis; + fthetady += Astate[aIndex + 2] * ibasis * jdbasis * kbasis; + fthetadz += Astate[aIndex + 2] * ibasis * jbasis * kdbasis; + break; } + } } + } } + + + if ((*configHash)["output_txt"] == "true"){ + samuraistream << std::scientific << i << "\t" << j << "\t" << k + << "\t" << pip << "\t" << thetarhop << "\t" << ftheta << "\t" << "\n"; + } + + // On the nodes + if (!ihalf and !jhalf and !khalf){ + int fIndex = (iDim-2)*(jDim-2)*(kDim-2); + int posIndex = (iDim-2)*(jDim-2)*(kIndex-1) + (iDim-2)*(jIndex-1) + (iIndex-1); + finalAnalysis[fIndex * 0 + posIndex] = pip; + finalAnalysis[fIndex * 1 + posIndex] = thetarhop; + finalAnalysis[fIndex * 2 + posIndex] = ftheta; + finalAnalysis[fIndex * 3 + posIndex] = pipdx; + finalAnalysis[fIndex * 4 + posIndex] = thetarhopdx; + finalAnalysis[fIndex * 5 + posIndex] = fthetadx; + finalAnalysis[fIndex * 6 + posIndex] = pipdy; + finalAnalysis[fIndex * 7 + posIndex] = thetarhopdy; + finalAnalysis[fIndex * 8 + posIndex] = fthetady; + finalAnalysis[fIndex * 9 + posIndex] = pipdz; + finalAnalysis[fIndex * 10 + posIndex] = thetarhopdz; + finalAnalysis[fIndex * 11 + posIndex] = fthetadz; + } + } } + } } + } } + } } + } - QString fileName = "samurai_XYZ_" + samuraiMode + "_" + suffix; - QString outFileName = outputPath.absoluteFilePath(fileName); - - // Write the Obs to a summary text file - if (configHash->value("output_qc") == "true") { - QString qcout = "samurai_QC_" + samuraiMode + "_" + suffix + ".out"; - QString qcFileName = outputPath.absoluteFilePath(qcout); - ofstream qcstream(qcFileName.toAscii().data()); - ostream_iterator os(qcstream, "\t "); - *os++ = "Observation"; - *os++ = "Inverse Error"; - *os++ = "X"; - *os++ = "Y"; - *os++ = "Z"; - *os++ = "Type"; - *os++ = "Time"; - *os++ = "pip"; - *os++ = "thetarhop"; - *os++ = "ftheta"; - *os++ = "Analysis"; - *os++ = "Background"; - qcstream << endl; - qcstream.precision(10); - - ostream_iterator od(qcstream, "\t "); - for (int m = 0; m < mObs; m++) { - int mi = m*(obMetaSize+varDim*derivDim); - real i = obsVector[mi+2]; - real j = obsVector[mi+3]; - real k = obsVector[mi+4]; - real tempsum = 0; - int ii = (int)((i - iMin)*DIrecip); - int jj = (int)((j - jMin)*DJrecip); - int kk = (int)((k - kMin)*DKrecip); - real ibasis = 0; - real jbasis = 0; - real kbasis = 0; - for (int var = 0; var < varDim; var++) { - for (int d = 0; d < derivDim; d++) { - int wgt_index = mi + obMetaSize +d*varDim + var; - if (!obsVector[wgt_index]) continue; - for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { - kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, derivative[d][2], kBCL[var], kBCR[var]); - for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { - int iNode = iiNode; - if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; - if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); - if ((iNode < 0) or (iNode >= iDim)) continue; - ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, derivative[d][1], iBCL[var], iBCR[var]); - - for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { - int jNode = jjNode; - if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; - if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); - if ((jNode < 0) or (jNode >= jDim)) continue; - int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; - jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, derivative[d][0], jBCL[var], jBCR[var]); - tempsum += Astate[aIndex + var] * ibasis * jbasis * kbasis * obsVector[wgt_index]; - } - } - } - } - } - for (int t=0; t < obMetaSize-1; t++) { - *od++ = obsVector[mi+t]; - } - int unixtime = (int)obsVector[mi+6]; - QDateTime obtime; - obtime.setTime_t(unixtime); - obtime.setTimeSpec(Qt::UTC); - QString timestring = obtime.toString("hh:mm:ss.zzz"); - qcstream << timestring.toStdString() << "\t"; - - // Multiply the weight by the ob -- Observations.in has individual weights already - // Only non-derivative for now - for (int t=obMetaSize; t os(qcstream, "\t "); + *os++ = "Observation"; + *os++ = "Inverse Error"; + *os++ = "X"; + *os++ = "Y"; + *os++ = "Z"; + *os++ = "Type"; + *os++ = "Time"; + *os++ = "pip"; + *os++ = "thetarhop"; + *os++ = "ftheta"; + *os++ = "Analysis"; + *os++ = "Background"; + qcstream << std::endl; + qcstream.precision(10); + std::ostream_iterator od(qcstream, "\t "); + for (int m = 0; m < mObs; m++) { + int mi = m*(obMetaSize+varDim*derivDim); + real i = obsVector[mi+2]; + real j = obsVector[mi+3]; + real k = obsVector[mi+4]; + real tempsum = 0; + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0; + real jbasis = 0; + real kbasis = 0; + for (int var = 0; var < varDim; var++) { + for (int d = 0; d < derivDim; d++) { + int wgt_index = mi + obMetaSize +d*varDim + var; + if (!obsVector[wgt_index]) continue; + for (int kNode = std::max(kk-1,0); kNode <= std::min(kk+2,kDim-1); ++kNode) { + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, derivative[d][2], kBCL[var], kBCR[var]); + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, derivative[d][1], iBCL[var], iBCR[var]); + + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, derivative[d][0], jBCL[var], jBCR[var]); + tempsum += Astate[aIndex + var] * ibasis * jbasis * kbasis * obsVector[wgt_index]; + } + } + } } - } + } + for (int t=0; t < obMetaSize-1; t++) { + *od++ = obsVector[mi+t]; + } + int unixtime = (int)obsVector[mi+6]; + std::time_t obtime = static_cast(unixtime); + std::tm* timeinfo = std::gmtime(&obtime); + char timestring[100]; + std::strftime(timestring, sizeof(timestring), "%H:%M:%S", timeinfo); + qcstream << timestring << "\t"; - adjustInternalDomain(-1); + // Multiply the weight by the ob -- Observations.in has individual weights already + // Only non-derivative for now + for (int t=obMetaSize; tvalue("output_netcdf") == "true") { - QString cdfFileName = outFileName + ".nc"; - if (!writeNetCDF(outputPath.absoluteFilePath(cdfFileName))) - cout << "Error writing netcdf file " << cdfFileName.toStdString() << endl; - } - // Write out to an asi file - if (configHash->value("output_asi") == "true") { - QString asiFileName = outFileName + ".asi"; - if (!writeAsi(outputPath.absoluteFilePath(asiFileName))) - cout << "Error writing asi file " << asiFileName.toStdString() << endl; } - // Set the domain back - adjustInternalDomain(1); + } - // Free the memory for the analysis variables - delete[] finalAnalysis; + adjustInternalDomain(-1); - return true; + // Write out to a netCDF file + if ((*configHash)["output_netcdf"] == "true") { + std::string cdfFileName = outFileName + ".nc"; + if (!writeNetCDF(outputPath + cdfFileName)) + std::cout << "Error writing netcdf file " << cdfFileName << std::endl; + } + // Write out to an asi file + if ((*configHash)["output_asi"] == "true") { + std::string asiFileName = outFileName + ".asi"; + if (!writeAsi(outputPath + asiFileName)) + std::cout << "Error writing asi file " << asiFileName << std::endl; + } + // Set the domain back + adjustInternalDomain(1); + // Free the memory for the analysis variables + delete[] finalAnalysis; + + return true; } bool CostFunctionXYZ::writeNetCDF(const std::string& netcdfFileName) From cd808b2c49cef287150a88bb5c6ea3ec7b4c66aa Mon Sep 17 00:00:00 2001 From: cenamiller Date: Tue, 20 Aug 2024 15:09:53 -0600 Subject: [PATCH 3/3] Removing readNetCDF_thermo(), calc_A_thermo()-calc_E_thermo(), getValue_thermo(), and getDerivative_thermo() from VarDriver. Added NetCDF_XYZ class from thermo code instead. Removed QTStrings and updated to NC3 --- src/CMakeLists.txt | 6 +- src/NetCDF.cpp | 33 ++++ src/NetCDF.h | 46 +++++ src/NetCDF_XYZ.cpp | 450 +++++++++++++++++++++++++++++++++++++++++++++ src/NetCDF_XYZ.h | 63 +++++++ src/VarDriver.cpp | 423 +----------------------------------------- src/VarDriver.h | 34 +--- 7 files changed, 603 insertions(+), 452 deletions(-) create mode 100644 src/NetCDF.cpp create mode 100644 src/NetCDF.h create mode 100644 src/NetCDF_XYZ.cpp create mode 100644 src/NetCDF_XYZ.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9d52d72..46c49c1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -80,7 +80,9 @@ set(samurai_HDRS Dorade.h ErrorData.h FrameCenter.h - MetObs.h + MetObs.h + NetCDF.h + NetCDF_XYZ.h Observation.h precision.h Projection.h @@ -116,6 +118,8 @@ set(common_SRCS MetObs.cpp mac_debug.xcconfig mac_release.xcconfig + NetCDF.cpp + NetCDF_XYZ.cpp Observation.cpp Projection.cpp RecursiveFilter.cpp diff --git a/src/NetCDF.cpp b/src/NetCDF.cpp new file mode 100644 index 0000000..89f2545 --- /dev/null +++ b/src/NetCDF.cpp @@ -0,0 +1,33 @@ +/* + * NetCDF.cpp + * samurai + * + * Created by Annette Foerster on 8/15/13. + * Adapted from $Id: pres_temp_4D_rd.cpp,v 1.11 2006/08/22 19:22:06 ed Exp $ + * (http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial) + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + + +#include "NetCDF.h" +#include + +NetCDF::NetCDF(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim) :c_p(1005.7), g(9.81), f(5.85e-05), pi(3.14159265358979323846) + { + NDIMS = 4; + NALT = metFile_kdim; + NX = metFile_idim; // NRADIUS or NLON + NY = metFile_jdim; // NTHETA or NLAT + NREC = 0; + NC_ERR = 2; + + std::cout << "Constructor \n"; +} + +NetCDF::~NetCDF() { +} + + + + diff --git a/src/NetCDF.h b/src/NetCDF.h new file mode 100644 index 0000000..f638cd1 --- /dev/null +++ b/src/NetCDF.h @@ -0,0 +1,46 @@ + /* + * NetCDF.h + * samurai + * + * Created by Annette Foerster on 8/15/13. + * Based on example code from netCDF library + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + +#ifndef NETCDF_H +#define NETCDF_H + +//#include +#include +#include + +class NetCDF +{ + +public: + NetCDF(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim); + + virtual ~NetCDF(); + + virtual int readNetCDF(const char* filename)=0; + virtual double getValue(const int &i,const int &j,const int &k,const std::string& varName)=0; + virtual double getDerivative(const int &i,const int &j,const int &k, const std::string &var, const int &der)=0; + virtual double calc_A(const int &i,const int &j,const int &k)=0; + virtual double calc_B(const int &i,const int &j,const int &k)=0; + virtual double calc_C(const int &i,const int &j,const int &k)=0; + virtual double calc_D(const int &i,const int &j,const int &k)=0; + virtual double calc_E(const int &i,const int &j,const int &k)=0; + + +protected: + int NDIMS, NALT, NX, NY, NREC, NC_ERR; + + const float c_p; + const float g; + const double f; //this needs to be made dynamic eventually, here lat assumed to be 22 deg north + const double pi; +}; + + +#endif diff --git a/src/NetCDF_XYZ.cpp b/src/NetCDF_XYZ.cpp new file mode 100644 index 0000000..eff84eb --- /dev/null +++ b/src/NetCDF_XYZ.cpp @@ -0,0 +1,450 @@ +/* + * NetCDF_RTZ.cpp + * samurai + * + * Created by Annette Foerster. + * Adapted from $Id: pres_temp_4D_rd.cpp,v 1.11 2006/08/22 19:22:06 ed Exp $ + * (http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial) + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + +#include "NetCDF.h" +#include "NetCDF_XYZ.h" +#include +#include +//#include +#include + + +NetCDF_XYZ::NetCDF_XYZ(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim) +:NetCDF::NetCDF(metFile_idim, metFile_jdim, metFile_kdim) + { + NLON = metFile_idim; + NLAT = metFile_jdim; + NALT = metFile_kdim; + + longitude = new float [NLON]; + latitude = new float [NLAT]; + altitude = new float [NALT]; + u = new float[NALT*NLON*NLAT]; + v = new float[NALT*NLON*NLAT]; + w = new float[NALT*NLON*NLAT]; + dudx = new float[NALT*NLON*NLAT]; + dvdx = new float[NALT*NLON*NLAT]; + dwdx = new float[NALT*NLON*NLAT]; + dudy = new float[NALT*NLON*NLAT]; + dvdy = new float[NALT*NLON*NLAT]; + dwdy = new float[NALT*NLON*NLAT]; + dudz = new float[NALT*NLON*NLAT]; + dvdz = new float[NALT*NLON*NLAT]; + dwdz = new float[NALT*NLON*NLAT]; + thetarhobar = new float[NALT*NLON*NLAT]; + dpibardx = new float[NALT*NLON*NLAT]; + dpibardy = new float[NALT*NLON*NLAT]; + + } + +NetCDF_XYZ::~NetCDF_XYZ() +{ + delete[] longitude; + delete[] latitude; + delete[] altitude; + delete[] u; + delete[] v; + delete[] w; + delete[] dudx; + delete[] dvdx; + delete[] dwdx; + delete[] dudy; + delete[] dvdy; + delete[] dwdy; + delete[] dudz; + delete[] dvdz; + delete[] dwdz; + delete[] thetarhobar; + delete[] dpibardx; + delete[] dpibardy; + +} + +int NetCDF_XYZ::readNetCDF(const char* filename) { + Nc3Error err(Nc3Error::verbose_nonfatal); + Nc3File dataFile(filename, Nc3File::ReadOnly); + + if(!dataFile.is_valid()) + return NC_ERR; + + // Get pointers to the latitude and longitude variables. + Nc3Var *lonVar, *latVar, *altVar; + if (!(lonVar = dataFile.get_var("lon"))) + return NC_ERR; + if (!(latVar = dataFile.get_var("lat"))) + return NC_ERR; + if (!(altVar = dataFile.get_var("z"))) + return NC_ERR; + + // Get the lat/lon data from the file. + if (!lonVar->get(longitude, NLON)) + return NC_ERR; + if (!latVar->get(latitude, NLAT)) + return NC_ERR; + if (!altVar->get(altitude, NALT)) + return NC_ERR; + + // Get pointers to the pressure and temperature variables. + Nc3Var *uVar,*vVar,*wVar,*dudxVar,*dvdxVar,*dwdxVar,*dudyVar,*dvdyVar,*dwdyVar,*dudzVar,*dvdzVar,*dwdzVar, *trbVar, *dpibdxVar, *dpibdyVar; + + if (!(uVar = dataFile.get_var("u"))) + return NC_ERR; + if (!(vVar = dataFile.get_var("v"))) + return NC_ERR; + if (!(wVar = dataFile.get_var("w"))) + return NC_ERR; + if (!(dudxVar = dataFile.get_var("dudx"))) + return NC_ERR; + if (!(dvdxVar = dataFile.get_var("dvdx"))) + return NC_ERR; + if (!(dwdxVar = dataFile.get_var("dwdx"))) + return NC_ERR; + if (!(dudyVar = dataFile.get_var("dudy"))) + return NC_ERR; + if (!(dvdyVar = dataFile.get_var("dvdy"))) + return NC_ERR; + if (!(dwdyVar = dataFile.get_var("dwdy"))) + return NC_ERR; + if (!(dudzVar = dataFile.get_var("dudz"))) + return NC_ERR; + if (!(dvdzVar = dataFile.get_var("dvdz"))) + return NC_ERR; + if (!(dwdzVar = dataFile.get_var("dwdz"))) + return NC_ERR; + if (!(trbVar = dataFile.get_var("trb"))) + return NC_ERR; + if (!(dpibdxVar = dataFile.get_var("dpibdx"))) + return NC_ERR; + if (!(dpibdyVar = dataFile.get_var("dpibdy"))) + return NC_ERR; + + if (!uVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!vVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!wVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dudxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dvdxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dwdxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dudyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dvdyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dwdyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dudzVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dvdzVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dwdzVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!trbVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dpibdxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dpibdyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + + if (!uVar->get(u, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!vVar->get(v, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!wVar->get(w, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dudxVar->get(dudx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dvdxVar->get(dvdx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dwdxVar->get(dwdx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dudyVar->get(dudy, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dvdyVar->get(dvdy, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dwdyVar->get(dwdy, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dudzVar->get(dudz, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dvdzVar->get(dvdz, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dwdzVar->get(dwdz, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!trbVar->get(thetarhobar, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dpibdxVar->get(dpibardx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dpibdyVar->get(dpibardy, 1, NALT, NLAT, NLON)) + return NC_ERR; + + return 0; +} + + +double NetCDF_XYZ::getValue(const int &i,const int &j,const int &k, const std::string& varName) +{ + //Returned values are all in SI units + double value_out; + + if ((i<0) or (i > NLON-1)) return false; + if ((j<0) or (j > NLAT-1)) return false; + if ((k<0) or (k > NALT-1)) return false; + + if (varName=="lon") { + value_out= longitude[i]; + } else if (varName=="lat") { + value_out= latitude[j]; + } else if (varName=="z") { + value_out= altitude[k]*1000.0; + } else if (varName=="u") { + value_out= u[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="v") { + value_out= v[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="w") { + value_out= w[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="dudx") { + value_out= dudx[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dvdx") { + value_out= dvdx[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dwdx") { + value_out= dwdx[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dudy") { + value_out= dudy[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dvdy") { + value_out= dvdy[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dwdy") { + value_out= dwdy[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dudz") { + value_out= dudz[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dvdz") { + value_out= dvdz[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dwdz") { + value_out= dwdz[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="trb") { + value_out= thetarhobar[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="dpibdx") { + value_out= dpibardx[k*NLAT*NLON + j*NLAT + i]/1000.0; + } else if (varName=="dpibdy") { + value_out= dpibardy[k*NLAT*NLON + j*NLAT + i]/1000.0; + } else if (varName=="A") { + value_out= this->calc_A(i,j,k); + } else if (varName=="B") { + value_out= this->calc_B(i,j,k); + } else if (varName=="C") { + value_out= this->calc_C(i,j,k); + + } else { + std::cout << "Requested Variable unknown. \n"; + std::cout << varName << "\n"; + return false; + } + + return value_out; +} + + +double NetCDF_XYZ::calc_A(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,"trb"); + double u = this->getValue(i,j,k,"u"); + double dudx = this->getValue(i,j,k,"dudx"); + double v = this->getValue(i,j,k,"v"); + double dudy = this->getValue(i,j,k,"dudy"); + double w = this->getValue(i,j,k,"w"); + double dudz = this->getValue(i,j,k,"dudz"); + double dpibdx = this->getValue(i,j,k,"dpibdx"); + float c_p = 1005.7; + + if (thetarhobar==-999 or u==-999 or dudx*1.0E5==-999 or v==-999 or dudy*1.0E5==-999 or w==-999 or dudz*1.0E5==-999 or dpibdx*1000.0==-999){ + return -999;} + + double a = 1.0/(c_p*thetarhobar)* (u*dudx+v*dudy+w*dudz-f*v)+dpibdx; + return a; +} + + +double NetCDF_XYZ::calc_B(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,"trb"); + double u = this->getValue(i,j,k,"u"); + double dvdx = this->getValue(i,j,k,"dvdx"); + double v = this->getValue(i,j,k,"v"); + double dvdy = this->getValue(i,j,k,"dvdy"); + double w = this->getValue(i,j,k,"w"); + double dvdz = this->getValue(i,j,k,"dvdz"); + double dpibdy = this->getValue(i,j,k,"dpibdy"); + float c_p = 1005.7; + + if (thetarhobar==-999 or u==-999 or dvdx*1.0E5==-999 or v==-999 or dvdy*1.0E5==-999 or w==-999 or dvdz*1.0E5==-999 or dpibdy*1000.0==-999){ + return -999;} + + double b = 1.0/(c_p*thetarhobar)*(u*dvdx+v*dvdy+w*dvdz+f*u)+dpibdy; // trp neglected + return b; +} + +double NetCDF_XYZ::calc_C(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,"trb"); + double u = this->getValue(i,j,k,"u"); + double dwdx = this->getValue(i,j,k,"dwdx"); + double v = this->getValue(i,j,k,"v"); + double dwdy = this->getValue(i,j,k,"dwdy"); + double w = this->getValue(i,j,k,"w"); + double dwdz = this->getValue(i,j,k,"dwdz"); + float c_p = 1005.7; + float g = 9.81; + + if (thetarhobar==-999 or u==-999 or dwdx*1.0E5==-999 or v==-999 or dwdy*1.0E5==-999 or w==-999 or dwdz*1.0E5==-999){ + return -999;} + + double c = 1.0/(c_p*thetarhobar)*(u*dwdx+v*dwdy+w*dwdz); + return c; +} + +double NetCDF_XYZ::calc_D(const int &i,const int &j,const int &k) +{ + + double thetarhobar = this->getValue(i,j,k,"trb"); + double dAdz = this->getDerivative(i,j,k,"A",3)*1000.0; // per km + double dCdx = this->getDerivative(i,j,k,"C",1)*1000.0; // per km + float c_p = 1005.7; + float g = 9.81; + + if (thetarhobar==-999 or dAdz==-999000 or dCdx==-999000){ + return -999;} + + double d = (dAdz-dCdx)*thetarhobar*thetarhobar*(-c_p/g); + + return d; +} + +double NetCDF_XYZ::calc_E(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,"trb"); + double dBdz = this->getDerivative(i,j,k,"B",3)*1000.0; // per km + double dCdy = this->getDerivative(i,j,k,"C",2)*1000.0; // per km + float c_p = 1005.7; + float g = 9.81; + + if (thetarhobar==-999 or dBdz==-999000 or dCdy==-999000){ + return -999;} + + double e = (dBdz-dCdy)*thetarhobar*thetarhobar*(-c_p/g); + + return e; +} + + +double NetCDF_XYZ::getDerivative(const int &i,const int &j,const int &k, const std::string &var, const int &der) +{ + // input variable "der" specifies direction of derivation: 1=dx, 2=dy, 3=dz + double derivative; + std::string derDir; + // Geographic functions + GeographicLib::TransverseMercatorExact tm = GeographicLib::TransverseMercatorExact::UTM(); + double referenceLon = -90.0; //arbitrary + double x1,x2,y1,y2; + + switch ( der ) { + case 1: + derDir = "X"; + if (i==0){ + //need to convert lat/lon differences to kms + tm.Forward(referenceLon,this->getValue(i+1,j,k,"lat"),this->getValue(i+1,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i+1,j,k,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i+1,j,k,var)-this->getValue(i,j,k,var))/distance; + } + } else if (i== NLON-1) { + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i-1,j,k,"lat"),this->getValue(i-1,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j,k,var)==-999 || this->getValue(i-1,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k,var)-this->getValue(i-1,j,k,var))/distance; + } + } else { + tm.Forward(referenceLon,this->getValue(i+1,j,k,"lat"),this->getValue(i+1,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i-1,j,k,"lat"),this->getValue(i-1,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i+1,j,k,var)==-999 || this->getValue(i-1,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i+1,j,k,var)-this->getValue(i-1,j,k,var))/distance; + } + } + break; + case 2: + derDir = "Y"; + if (j==0){ + tm.Forward(referenceLon,this->getValue(i,j+1,k,"lat"),this->getValue(i,j+1,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j+1,k,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j+1,k,var)-this->getValue(i,j,k,var))/distance; + } + } else if (j== NLAT-1) { + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i,j-1,k,"lat"),this->getValue(i,j-1,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j,k,var)==-999 || this->getValue(i,j-1,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k,var)-this->getValue(i,j-1,k,var))/distance; + } + } else { + tm.Forward(referenceLon,this->getValue(i,j+1,k,"lat"),this->getValue(i,j+1,k,"lon"),x2,y2); + tm.Forward(referenceLon,this->getValue(i,j-1,k,"lat"),this->getValue(i,j-1,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j+1,k,var)==-999 || this->getValue(i,j-1,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j+1,k,var)-this->getValue(i,j-1,k,var))/distance; + } + } + break; + case 3: + derDir = "Z"; + if (k==0){ + if (this->getValue(i,j,k+1,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k+1,var)-this->getValue(i,j,k,var))/(this->getValue(i,j,k+1,"z")-this->getValue(i,j,k,"z")); + } + } else if (k== NALT-1) { + if (this->getValue(i+1,j,k,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k,var)-this->getValue(i,j,k-1,var))/(this->getValue(i,j,k,"z")-this->getValue(i,j,k-1,"z")); + } + } else { + if (this->getValue(i,j,k+1,var)==-999 || this->getValue(i,j,k-1,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k+1,var)-this->getValue(i,j,k-1,var))/(this->getValue(i,j,k+1,"z")-this->getValue(i,j,k-1,"z")); + } + } + break; + default: + std::cout << "Unknown value for calculating derivative. Valid options are 1, 2 and 3.\n"; + exit(1); + } + return derivative; +} diff --git a/src/NetCDF_XYZ.h b/src/NetCDF_XYZ.h new file mode 100644 index 0000000..2df6247 --- /dev/null +++ b/src/NetCDF_XYZ.h @@ -0,0 +1,63 @@ + /* + * NetCDF.h + * samurai + * + * Created by Annette Foerster. + * Based on example code from netCDF library + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + +#ifndef NETCDF_XYZ_H +#define NETCDF_XYZ_H + +#include "NetCDF.h" +//#include +#include +#include + +class NetCDF_XYZ : public NetCDF +{ + +public: + NetCDF_XYZ(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim); + + ~NetCDF_XYZ(); + + int readNetCDF(const char* filename); + double getValue(const int &i,const int &j,const int &k,const std::string& varName); + double getDerivative(const int &i,const int &j,const int &k, const std::string &var, const int &der); + double calc_A(const int &i,const int &j,const int &k); + double calc_B(const int &i,const int &j,const int &k); + double calc_C(const int &i,const int &j,const int &k); + double calc_D(const int &i,const int &j,const int &k); + double calc_E(const int &i,const int &j,const int &k); + +protected: + int NLON, NLAT; + + float* longitude; + float* latitude; + float* altitude; + float* u; + float* v; + float* w; + float* dudx; + float* dvdx; + float* dwdx; + float* dudy; + float* dvdy; + float* dwdy; + float* dudz; + float* dvdz; + float* dwdz; + float* thetarhobar; + float* dpibardx; + float* dpibardy; + + std::string varName; + +}; + + +#endif diff --git a/src/VarDriver.cpp b/src/VarDriver.cpp index 3204803..1b311f6 100644 --- a/src/VarDriver.cpp +++ b/src/VarDriver.cpp @@ -26,7 +26,6 @@ #include // Constructor VarDriver::VarDriver() - :f_thermo(5.85e-05) // this needs to be made dynamic eventually, here lat assumed to be 22 deg north { // Constant for all drivers Pi = acos(-1); @@ -64,53 +63,11 @@ VarDriver::VarDriver() // By default we have fixed grid dimensions coming from the config file fixedGrid = true; - // enable for thermo - if(0){ - longitude = new float [NLON]; - latitude = new float [NLAT]; - altitude = new float [NALT]; - u = new float[NALT*NLON*NLAT]; - v = new float[NALT*NLON*NLAT]; - w = new float[NALT*NLON*NLAT]; - dudx = new float[NALT*NLON*NLAT]; - dvdx = new float[NALT*NLON*NLAT]; - dwdx = new float[NALT*NLON*NLAT]; - dudy = new float[NALT*NLON*NLAT]; - dvdy = new float[NALT*NLON*NLAT]; - dwdy = new float[NALT*NLON*NLAT]; - dudz = new float[NALT*NLON*NLAT]; - dvdz = new float[NALT*NLON*NLAT]; - dwdz = new float[NALT*NLON*NLAT]; - thetarhobar = new float[NALT*NLON*NLAT]; - dpibardx = new float[NALT*NLON*NLAT]; - dpibardy = new float[NALT*NLON*NLAT]; - } } // Destructor VarDriver::~VarDriver() { - //enable for thermo - if(0){ - delete[] longitude; - delete[] latitude; - delete[] altitude; - delete[] u; - delete[] v; - delete[] w; - delete[] dudx; - delete[] dvdx; - delete[] dwdx; - delete[] dudy; - delete[] dvdy; - delete[] dwdy; - delete[] dudz; - delete[] dvdz; - delete[] dwdz; - delete[] thetarhobar; - delete[] dpibardx; - delete[] dpibardy; - } } @@ -1987,385 +1944,7 @@ Projection::ProjectionType VarDriver::projectionFromConfig() return retVal; } -// This routine from the thermo branch, NetCDF_XYZ::readNetCDF - -int VarDriver::readNetCDF_thermo(const char* filename){ - - - Nc3Error err(Nc3Error::verbose_nonfatal); - Nc3File dataFile(filename, Nc3File::ReadOnly); - - if(!dataFile.is_valid()) - return NC_ERR; - - // Get pointers to the latitude and longitude variables. - Nc3Var *lonVar, *latVar, *altVar; - if (!(lonVar = dataFile.get_var("lon"))) - return NC_ERR; - if (!(latVar = dataFile.get_var("lat"))) - return NC_ERR; - if (!(altVar = dataFile.get_var("z"))) - return NC_ERR; - - // Get the lat/lon data from the file. - if (!lonVar->get(longitude, NLON)) - return NC_ERR; - if (!latVar->get(latitude, NLAT)) - return NC_ERR; - if (!altVar->get(altitude, NALT)) - return NC_ERR; - - // Get pointers to the pressure and temperature variables. - Nc3Var *uVar,*vVar,*wVar,*dudxVar,*dvdxVar,*dwdxVar,*dudyVar,*dvdyVar,*dwdyVar,*dudzVar,*dvdzVar,*dwdzVar, *trbVar, *dpibdxVar, *dpibdyVar; - - if (!(uVar = dataFile.get_var("u"))) - return NC_ERR; - if (!(vVar = dataFile.get_var("v"))) - return NC_ERR; - if (!(wVar = dataFile.get_var("w"))) - return NC_ERR; - if (!(dudxVar = dataFile.get_var("dudx"))) - return NC_ERR; - if (!(dvdxVar = dataFile.get_var("dvdx"))) - return NC_ERR; - if (!(dwdxVar = dataFile.get_var("dwdx"))) - return NC_ERR; - if (!(dudyVar = dataFile.get_var("dudy"))) - return NC_ERR; - if (!(dvdyVar = dataFile.get_var("dvdy"))) - return NC_ERR; - if (!(dwdyVar = dataFile.get_var("dwdy"))) - return NC_ERR; - if (!(dudzVar = dataFile.get_var("dudz"))) - return NC_ERR; - if (!(dvdzVar = dataFile.get_var("dvdz"))) - return NC_ERR; - if (!(dwdzVar = dataFile.get_var("dwdz"))) - return NC_ERR; - if (!(trbVar = dataFile.get_var("trb"))) - return NC_ERR; - if (!(dpibdxVar = dataFile.get_var("dpibdx"))) - return NC_ERR; - if (!(dpibdyVar = dataFile.get_var("dpibdy"))) - return NC_ERR; - - if (!uVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!vVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!wVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dudxVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dvdxVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dwdxVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dudyVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dvdyVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dwdyVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dudzVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dvdzVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dwdzVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!trbVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dpibdxVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - if (!dpibdyVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; - - if (!uVar->get(u, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!vVar->get(v, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!wVar->get(w, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dudxVar->get(dudx, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dvdxVar->get(dvdx, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dwdxVar->get(dwdx, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dudyVar->get(dudy, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dvdyVar->get(dvdy, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dwdyVar->get(dwdy, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dudzVar->get(dudz, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dvdzVar->get(dvdz, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dwdzVar->get(dwdz, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!trbVar->get(thetarhobar, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dpibdxVar->get(dpibardx, 1, NALT, NLAT, NLON)) - return NC_ERR; - if (!dpibdyVar->get(dpibardy, 1, NALT, NLAT, NLON)) - return NC_ERR; - - return 0; -} - -double VarDriver::getValue_thermo(const int &i,const int &j,const int &k, const std::string& varName){ - //Returned values are all in SI units - double value_out; - - if ((i<0) or (i > NLON-1)) return false; - if ((j<0) or (j > NLAT-1)) return false; - if ((k<0) or (k > NALT-1)) return false; - - if (varName=="lon") { - value_out= longitude[i]; - } else if (varName=="lat") { - value_out= latitude[j]; - } else if (varName=="z") { - value_out= altitude[k]*1000.0; - } else if (varName=="u") { - value_out= u[k*NLAT*NLON + j*NLAT + i]; - } else if (varName=="v") { - value_out= v[k*NLAT*NLON + j*NLAT + i]; - } else if (varName=="w") { - value_out= w[k*NLAT*NLON + j*NLAT + i]; - } else if (varName=="dudx") { - value_out= dudx[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dvdx") { - value_out= dvdx[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dwdx") { - value_out= dwdx[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dudy") { - value_out= dudy[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dvdy") { - value_out= dvdy[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dwdy") { - value_out= dwdy[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dudz") { - value_out= dudz[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dvdz") { - value_out= dvdz[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="dwdz") { - value_out= dwdz[k*NLAT*NLON + j*NLAT + i]/1.0E5; - } else if (varName=="trb") { - value_out= thetarhobar[k*NLAT*NLON + j*NLAT + i]; - } else if (varName=="dpibdx") { - value_out= dpibardx[k*NLAT*NLON + j*NLAT + i]/1000.0; - } else if (varName=="dpibdy") { - value_out= dpibardy[k*NLAT*NLON + j*NLAT + i]/1000.0; - } else if (varName=="A") { - value_out= this->calc_A_thermo(i,j,k); - } else if (varName=="B") { - value_out= this->calc_B_thermo(i,j,k); - } else if (varName=="C") { - value_out= this->calc_C_thermo(i,j,k); - - } else { - std::cout << "Requested Variable unknown. \n"; - std::cout << varName << "\n"; - return false; - } - - return value_out; -} - -double VarDriver::calc_A_thermo(const int &i,const int &j,const int &k) -{ - double thetarhobar = this->getValue_thermo(i,j,k,"trb"); - double u = this->getValue_thermo(i,j,k,"u"); - double dudx = this->getValue_thermo(i,j,k,"dudx"); - double v = this->getValue_thermo(i,j,k,"v"); - double dudy = this->getValue_thermo(i,j,k,"dudy"); - double w = this->getValue_thermo(i,j,k,"w"); - double dudz = this->getValue_thermo(i,j,k,"dudz"); - double dpibdx = this->getValue_thermo(i,j,k,"dpibdx"); - float c_p = 1005.7; - - if (thetarhobar==-999 or u==-999 or dudx*1.0E5==-999 or v==-999 or dudy*1.0E5==-999 or w==-999 or dudz*1.0E5==-999 or dpibdx*1000.0==-999){ - return -999;} - - double a = 1.0/(c_p*thetarhobar)* (u*dudx+v*dudy+w*dudz-f_thermo*v)+dpibdx; - return a; -} - -double VarDriver::calc_B_thermo(const int &i,const int &j,const int &k) -{ - double thetarhobar = this->getValue_thermo(i,j,k,"trb"); - double u = this->getValue_thermo(i,j,k,"u"); - double dvdx = this->getValue_thermo(i,j,k,"dvdx"); - double v = this->getValue_thermo(i,j,k,"v"); - double dvdy = this->getValue_thermo(i,j,k,"dvdy"); - double w = this->getValue_thermo(i,j,k,"w"); - double dvdz = this->getValue_thermo(i,j,k,"dvdz"); - double dpibdy = this->getValue_thermo(i,j,k,"dpibdy"); - float c_p = 1005.7; - - if (thetarhobar==-999 or u==-999 or dvdx*1.0E5==-999 or v==-999 or dvdy*1.0E5==-999 or w==-999 or dvdz*1.0E5==-999 or dpibdy*1000.0==-999){ - return -999;} - - double b = 1.0/(c_p*thetarhobar)*(u*dvdx+v*dvdy+w*dvdz+f_thermo*u)+dpibdy; // trp neglected - return b; -} - -double VarDriver::calc_C_thermo(const int &i,const int &j,const int &k) -{ - double thetarhobar = this->getValue_thermo(i,j,k,"trb"); - double u = this->getValue_thermo(i,j,k,"u"); - double dwdx = this->getValue_thermo(i,j,k,"dwdx"); - double v = this->getValue_thermo(i,j,k,"v"); - double dwdy = this->getValue_thermo(i,j,k,"dwdy"); - double w = this->getValue_thermo(i,j,k,"w"); - double dwdz = this->getValue_thermo(i,j,k,"dwdz"); - float c_p = 1005.7; - float g = 9.81; - - if (thetarhobar==-999 or u==-999 or dwdx*1.0E5==-999 or v==-999 or dwdy*1.0E5==-999 or w==-999 or dwdz*1.0E5==-999){ - return -999;} - - double c = 1.0/(c_p*thetarhobar)*(u*dwdx+v*dwdy+w*dwdz); - return c; -} - -double VarDriver::calc_D_thermo(const int &i,const int &j,const int &k) -{ - - double thetarhobar = this->getValue_thermo(i,j,k,"trb"); - double dAdz = this->getDerivative_thermo(i,j,k,"A",3)*1000.0; // per km - double dCdx = this->getDerivative_thermo(i,j,k,"C",1)*1000.0; // per km - float c_p = 1005.7; - float g = 9.81; - - if (thetarhobar==-999 or dAdz==-999000 or dCdx==-999000){ - return -999;} - - double d = (dAdz-dCdx)*thetarhobar*thetarhobar*(-c_p/g); - - return d; -} - -double VarDriver::calc_E_thermo(const int &i,const int &j,const int &k) -{ - double thetarhobar = this->getValue_thermo(i,j,k,"trb"); - double dBdz = this->getDerivative_thermo(i,j,k,"B",3)*1000.0; // per km - double dCdy = this->getDerivative_thermo(i,j,k,"C",2)*1000.0; // per km - float c_p = 1005.7; - float g = 9.81; - - if (thetarhobar==-999 or dBdz==-999000 or dCdy==-999000){ - return -999;} - - double e = (dBdz-dCdy)*thetarhobar*thetarhobar*(-c_p/g); - - return e; -} - -double VarDriver::getDerivative_thermo(const int &i,const int &j,const int &k, const std::string &var, const int &der) -{ - // input variable "der" specifies direction of derivation: 1=dx, 2=dy, 3=dz - double derivative; - std::string derDir; - // Geographic functions - GeographicLib::TransverseMercatorExact tm = GeographicLib::TransverseMercatorExact::UTM(); - double referenceLon = -90.0; //arbitrary - double x1,x2,y1,y2; - - switch ( der ) { - case 1: - derDir = "X"; - if (i==0){ - //need to convert lat/lon differences to kms - tm.Forward(referenceLon,this->getValue_thermo(i+1,j,k,"lat"),this->getValue_thermo(i+1,j,k,"lon"),x1,y1); - tm.Forward(referenceLon,this->getValue_thermo(i,j,k,"lat"),this->getValue_thermo(i,j,k,"lon"),x2,y2); - double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); - if (this->getValue_thermo(i+1,j,k,var)==-999 || this->getValue_thermo(i,j,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i+1,j,k,var)-this->getValue_thermo(i,j,k,var))/distance; - } - } else if (i== NLON-1) { - tm.Forward(referenceLon,this->getValue_thermo(i,j,k,"lat"),this->getValue_thermo(i,j,k,"lon"),x1,y1); - tm.Forward(referenceLon,this->getValue_thermo(i-1,j,k,"lat"),this->getValue_thermo(i-1,j,k,"lon"),x2,y2); - double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); - if (this->getValue_thermo(i,j,k,var)==-999 || this->getValue_thermo(i-1,j,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j,k,var)-this->getValue_thermo(i-1,j,k,var))/distance; - } - } else { - tm.Forward(referenceLon,this->getValue_thermo(i+1,j,k,"lat"),this->getValue_thermo(i+1,j,k,"lon"),x1,y1); - tm.Forward(referenceLon,this->getValue_thermo(i-1,j,k,"lat"),this->getValue_thermo(i-1,j,k,"lon"),x2,y2); - double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); - if (this->getValue_thermo(i+1,j,k,var)==-999 || this->getValue_thermo(i-1,j,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i+1,j,k,var)-this->getValue_thermo(i-1,j,k,var))/distance; - } - } - break; - case 2: - derDir = "Y"; - if (j==0){ - tm.Forward(referenceLon,this->getValue_thermo(i,j+1,k,"lat"),this->getValue_thermo(i,j+1,k,"lon"),x1,y1); - tm.Forward(referenceLon,this->getValue_thermo(i,j,k,"lat"),this->getValue_thermo(i,j,k,"lon"),x2,y2); - double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); - if (this->getValue_thermo(i,j+1,k,var)==-999 || this->getValue_thermo(i,j,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j+1,k,var)-this->getValue_thermo(i,j,k,var))/distance; - } - } else if (j== NLAT-1) { - tm.Forward(referenceLon,this->getValue_thermo(i,j,k,"lat"),this->getValue_thermo(i,j,k,"lon"),x1,y1); - tm.Forward(referenceLon,this->getValue_thermo(i,j-1,k,"lat"),this->getValue_thermo(i,j-1,k,"lon"),x2,y2); - double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); - if (this->getValue_thermo(i,j,k,var)==-999 || this->getValue_thermo(i,j-1,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j,k,var)-this->getValue_thermo(i,j-1,k,var))/distance; - } - } else { - tm.Forward(referenceLon,this->getValue_thermo(i,j+1,k,"lat"),this->getValue_thermo(i,j+1,k,"lon"),x2,y2); - tm.Forward(referenceLon,this->getValue_thermo(i,j-1,k,"lat"),this->getValue_thermo(i,j-1,k,"lon"),x2,y2); - double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); - if (this->getValue_thermo(i,j+1,k,var)==-999 || this->getValue_thermo(i,j-1,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j+1,k,var)-this->getValue_thermo(i,j-1,k,var))/distance; - } - } - break; - case 3: - derDir = "Z"; - if (k==0){ - if (this->getValue_thermo(i,j,k+1,var)==-999 || this->getValue_thermo(i,j,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j,k+1,var)-this->getValue_thermo(i,j,k,var))/(this->getValue_thermo(i,j,k+1,"z")-this->getValue_thermo(i,j,k,"z")); - } - } else if (k== NALT-1) { - if (this->getValue_thermo(i+1,j,k,var)==-999 || this->getValue_thermo(i,j,k,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j,k,var)-this->getValue_thermo(i,j,k-1,var))/(this->getValue_thermo(i,j,k,"z")-this->getValue_thermo(i,j,k-1,"z")); - } - } else { - if (this->getValue_thermo(i,j,k+1,var)==-999 || this->getValue_thermo(i,j,k-1,var)==-999){ - derivative = -999; - } else { - derivative = (this->getValue_thermo(i,j,k+1,var)-this->getValue_thermo(i,j,k-1,var))/(this->getValue_thermo(i,j,k+1,"z")-this->getValue_thermo(i,j,k-1,"z")); - } - } - break; - default: - std::cout << "Unknown value for calculating derivative. Valid options are 1, 2 and 3.\n"; - exit(1); - } - return derivative; -} +// bool VarDriver::read_mesonet(std::string& filename, std::vector* metObVector) { diff --git a/src/VarDriver.h b/src/VarDriver.h index f024194..f31529d 100644 --- a/src/VarDriver.h +++ b/src/VarDriver.h @@ -79,11 +79,11 @@ class VarDriver // bool readTerrain(std::string &filename, std::vector* metData); HashMap *getConfigHash() { return &configHash; } - double calc_A_thermo(const int &i,const int &j,const int &k); - double calc_B_thermo(const int &i,const int &j,const int &k); - double calc_C_thermo(const int &i,const int &j,const int &k); - double calc_D_thermo(const int &i,const int &j,const int &k); - double calc_E_thermo(const int &i,const int &j,const int &k); + //double calc_A_thermo(const int &i,const int &j,const int &k); + //double calc_B_thermo(const int &i,const int &j,const int &k); + //double calc_C_thermo(const int &i,const int &j,const int &k); + //double calc_D_thermo(const int &i,const int &j,const int &k); + //double calc_E_thermo(const int &i,const int &j,const int &k); protected: bool fixedGrid; // Indicates if the grid dims come from the config file or run call @@ -91,27 +91,6 @@ class VarDriver real CoriolisF; real Pi; int NC_ERR, NLON, NLAT, NALT, NREC; // Error code for netCDF from thermo - // netcdf thermo variables - float* longitude; - float* latitude; - float* altitude; - float* u; - float* v; - float* w; - float* dudx; - float* dvdx; - float* dwdx; - float* dudy; - float* dvdy; - float* dwdy; - float* dudz; - float* dvdz; - float* dwdz; - float* thetarhobar; - float* dpibardx; - float* dpibardy; - const double f_thermo; //this needs to be made dynamic eventually, here lat assumed to be 22 deg north - // end netcdf thermo variables unsigned int numVars; unsigned int numHeights; unsigned int maxHeights; @@ -188,9 +167,6 @@ class VarDriver bool readFrameCenters(); bool parseXMLconfig(const XMLNode& config); bool parseSamuraiConfig(const samurai_config &config); - int readNetCDF_thermo(const char* filename); - double getValue_thermo(const int &i,const int &j,const int &k, const std::string& varName); - double getDerivative_thermo(const int &i,const int &j,const int &k, const std::string &var, const int &der); Projection::ProjectionType projectionFromConfig(); };