diff --git a/ncar_scripts/TDRP/beltramiTR.tdrp b/ncar_scripts/TDRP/beltramiTR.tdrp index fa5763c..9c8503e 100644 --- a/ncar_scripts/TDRP/beltramiTR.tdrp +++ b/ncar_scripts/TDRP/beltramiTR.tdrp @@ -27,6 +27,14 @@ debug_kd = 0; debug_kd_step = 0; +///////////// wind_file ////////////////////////// +// +// Needed by Samurai THERMO mode +// +// Type: string + +wind_file = "../ncar_scripts/input/beltrami_test_v2.nc"; + //====================================================================== // // BACKGROUND SECTION. @@ -981,6 +989,36 @@ mesonet_rhow_error = 1; mesonet_tempk_error = 1; +///////////// thermo_A_error ///////////////////// +// +// Type: float + +thermo_A_error = 0.01; + +///////////// thermo_B_error ///////////////////// +// +// Type: float + +thermo_B_error = 0.01; + +///////////// thermo_C_error ///////////////////// +// +// Type: float + +thermo_C_error = 0.01; + +///////////// thermo_D_error ///////////////////// +// +// Type: float + +thermo_D_error = 0.01; + +///////////// thermo_E_error ///////////////////// +// +// Type: float + +thermo_E_error = 0.01; + //====================================================================== // // XYP SPECIFIC SECTION. @@ -1172,6 +1210,27 @@ bg_rhow_error = { 20, 2 }; bg_tempk_error = { 20, 2 }; +///////////// bg_pip_error ////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_pip_error = { 10, 10 }; + +///////////// bg_thetarhop_error ////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_thetarhop_error = { 10, 10 }; + +///////////// bg_ftheta_error ////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_ftheta_error = { 10, 10 }; + ///////////// i_filter_length ///////////////////////// // // Type: float @@ -1245,106 +1304,106 @@ k_max_wavenumber = { -1, -1 }; // // Type: int -i_pip_bcL = -1; +i_pip_bcL = "R0"; ///////////// i_pip_bcR /////////////////////////////// // // Type: int -i_pip_bcR = -1; +i_pip_bcR = "R0"; ///////////// i_thetarhop_bcL /////////////////////////////// // // Type: int -i_thetarhop_bcL = -1; +i_thetarhop_bcL = "R0"; ///////////// i_thetarhop_bcR /////////////////////////////// // // Type: int -i_thetarhop_bcR = -1; +i_thetarhop_bcR = "R0"; ///////////// i_ftheta_bcL /////////////////////////////// // // Type: int -i_ftheta_bcL = -1; +i_ftheta_bcL = "R0"; ///////////// i_ftheta_bcR /////////////////////////////// // // Type: int -i_ftheta_bcR = -1; +i_ftheta_bcR = "R0"; ///////////// j_pip_bcL /////////////////////////////// // // Type: int -j_pip_bcL = -1; +j_pip_bcL = "R0"; ///////////// j_pip_bcR /////////////////////////////// // // Type: int -j_pip_bcR = -1; +j_pip_bcR = "R0"; ///////////// j_thetarhop_bcL /////////////////////////////// // // Type: int -j_thetarhop_bcL = -1; +j_thetarhop_bcL = "R0"; ///////////// j_thetarhop_bcR /////////////////////////////// // // Type: int -j_thetarhop_bcR = -1; +j_thetarhop_bcR = "R0"; ///////////// j_ftheta_bcL /////////////////////////////// // // Type: int -j_ftheta_bcL = -1; +j_ftheta_bcL = "R0"; ///////////// j_ftheta_bcR /////////////////////////////// // // Type: int -j_ftheta_bcR = -1; +j_ftheta_bcR = "R0"; ///////////// k_pip_bcL /////////////////////////////// // // Type: int -k_pip_bcL = -1; +k_pip_bcL = "R0"; ///////////// k_pip_bcR /////////////////////////////// // // Type: int -k_pip_bcR = -1; +k_pip_bcR = "R0"; ///////////// k_thetarhop_bcL /////////////////////////////// // // Type: int -k_thetarhop_bcL = -1; +k_thetarhop_bcL = "R0"; ///////////// k_thetarhop_bcR /////////////////////////////// // // Type: int -k_thetarhop_bcR = -1; +k_thetarhop_bcR = "R0"; ///////////// k_ftheta_bcL /////////////////////////////// // // Type: int -k_ftheta_bcL = -1; +k_ftheta_bcL = "R0"; ///////////// k_ftheta_bcR /////////////////////////////// // // Type: int -k_ftheta_bcR = -1; \ No newline at end of file +k_ftheta_bcR = "R0"; diff --git a/ncar_scripts/input/beltrami_test_v2.nc b/ncar_scripts/input/beltrami_test_v2.nc new file mode 100644 index 0000000..6c29ead Binary files /dev/null and b/ncar_scripts/input/beltrami_test_v2.nc differ diff --git a/src/Args.cpp b/src/Args.cpp index c73e876..a3af808 100644 --- a/src/Args.cpp +++ b/src/Args.cpp @@ -101,6 +101,7 @@ bool Args::paramsToHash(HashMap *configHash) { CONFIG_INSERT_STR(debug_bgu_nc); CONFIG_INSERT_STR(debug_bgu_overwrite); CONFIG_INSERT_STR(fractl_nc_file); + CONFIG_INSERT_STR(wind_file); CONFIG_INSERT_BOOL(horizontal_radar_appx); CONFIG_INSERT_BOOL(load_background); CONFIG_INSERT_BOOL(load_bg_coefficients); @@ -280,13 +281,20 @@ bool Args::paramsToHash(HashMap *configHash) { CONFIG_INSERT_FLOAT(output_pressure_increment); CONFIG_INSERT_FLOAT(qscat_rhou_error); CONFIG_INSERT_FLOAT(qscat_rhov_error); + CONFIG_INSERT_FLOAT(thermo_A_error); + CONFIG_INSERT_FLOAT(thermo_B_error); + CONFIG_INSERT_FLOAT(thermo_C_error); + CONFIG_INSERT_FLOAT(thermo_D_error); + CONFIG_INSERT_FLOAT(thermo_E_error); CONFIG_INSERT_FLOAT(radar_fallspeed_error); CONFIG_INSERT_FLOAT(radar_min_error); CONFIG_INSERT_FLOAT(radar_sw_error); CONFIG_INSERT_FLOAT(rain_dbz); CONFIG_INSERT_FLOAT(sfmr_windspeed_error); + for (int iter = 1; iter <= params.num_iterations; iter++) { + // for WIND analysis CONFIG_INSERT_FLOAT_ARRAY(bg_qr_error, iter); CONFIG_INSERT_FLOAT_ARRAY(bg_rhoa_error, iter); CONFIG_INSERT_FLOAT_ARRAY(bg_rhou_error, iter); @@ -295,6 +303,11 @@ bool Args::paramsToHash(HashMap *configHash) { CONFIG_INSERT_FLOAT_ARRAY(bg_tempk_error, iter); CONFIG_INSERT_FLOAT_ARRAY(bg_qv_error, iter); + // for THERMO analysis + CONFIG_INSERT_FLOAT_ARRAY(bg_pip_error, iter); + CONFIG_INSERT_FLOAT_ARRAY(bg_thetarhop_error, iter); + CONFIG_INSERT_FLOAT_ARRAY(bg_ftheta_error, iter); + CONFIG_INSERT_FLOAT_ARRAY(i_filter_length, iter); CONFIG_INSERT_FLOAT_ARRAY(j_filter_length, iter); CONFIG_INSERT_FLOAT_ARRAY(k_filter_length, iter); @@ -309,26 +322,26 @@ bool Args::paramsToHash(HashMap *configHash) { CONFIG_INSERT_FLOAT_ARRAY(dirichlet_w_weight, iter); } // arguments required by THERMO - CONFIG_INSERT_INT(i_pip_bcL); - CONFIG_INSERT_INT(i_pip_bcR); - CONFIG_INSERT_INT(i_thetarhop_bcL); - CONFIG_INSERT_INT(i_thetarhop_bcR); - CONFIG_INSERT_INT(i_ftheta_bcL); - CONFIG_INSERT_INT(i_ftheta_bcR); - - CONFIG_INSERT_INT(j_pip_bcL); - CONFIG_INSERT_INT(j_pip_bcR); - CONFIG_INSERT_INT(j_thetarhop_bcL); - CONFIG_INSERT_INT(j_thetarhop_bcR); - CONFIG_INSERT_INT(j_ftheta_bcL); - CONFIG_INSERT_INT(j_ftheta_bcR); - - CONFIG_INSERT_INT(k_pip_bcL); - CONFIG_INSERT_INT(k_pip_bcR); - CONFIG_INSERT_INT(k_thetarhop_bcL); - CONFIG_INSERT_INT(k_thetarhop_bcR); - CONFIG_INSERT_INT(k_ftheta_bcL); - CONFIG_INSERT_INT(k_ftheta_bcR); + CONFIG_INSERT_STR(i_pip_bcL); + CONFIG_INSERT_STR(i_pip_bcR); + CONFIG_INSERT_STR(i_thetarhop_bcL); + CONFIG_INSERT_STR(i_thetarhop_bcR); + CONFIG_INSERT_STR(i_ftheta_bcL); + CONFIG_INSERT_STR(i_ftheta_bcR); + + CONFIG_INSERT_STR(j_pip_bcL); + CONFIG_INSERT_STR(j_pip_bcR); + CONFIG_INSERT_STR(j_thetarhop_bcL); + CONFIG_INSERT_STR(j_thetarhop_bcR); + CONFIG_INSERT_STR(j_ftheta_bcL); + CONFIG_INSERT_STR(j_ftheta_bcR); + + CONFIG_INSERT_STR(k_pip_bcL); + CONFIG_INSERT_STR(k_pip_bcR); + CONFIG_INSERT_STR(k_thetarhop_bcL); + CONFIG_INSERT_STR(k_thetarhop_bcR); + CONFIG_INSERT_STR(k_ftheta_bcL); + CONFIG_INSERT_STR(k_ftheta_bcR); return true; } diff --git a/src/CostFunction.cpp b/src/CostFunction.cpp index c049de1..8a2c31e 100644 --- a/src/CostFunction.cpp +++ b/src/CostFunction.cpp @@ -105,6 +105,7 @@ void CostFunction::truncatedNewton(real* qstate, real* g, const real ftol) real n_init_grad, n_grad, grad_dot; real gg; real rr, pAp, rr_m1, r_norm, r0_norm, rel_resid; + real r0_init_norm; real beta, alpha; real initstep; @@ -150,6 +151,7 @@ void CostFunction::truncatedNewton(real* qstate, real* g, const real ftol) if (its == 0) { f_init = f_val; n_init_grad = n_grad; + r0_init_norm = n_init_grad; //cout << "Newton: Initial Norm of the Gradient = " << n_init_grad << endl; } @@ -189,12 +191,14 @@ void CostFunction::truncatedNewton(real* qstate, real* g, const real ftol) r[j] = -g[j]; //initial residual r0 = b -Ax = -g p[j] = r[j]; //inital direction po = r0 } + //JMD: I am wondering if r0_norm should be the initial residual n_init_grad + // r0_norm = n_grad; r0_norm = n_grad; neg_curve = 0; //negative curvature check if (verbose) { - std::cout << "\t\tCG iteration " << std::setw(7) << std::right << "0" << ": r_norm = " << std::fixed << std::setw(20) << std::setprecision(8) << std::right << r0_norm << " rel_resid = " << std::setw(14) << std::setprecision(10) << std::right << 1.0 << endl; + std::cout << "\t\tCG iteration " << std::setw(7) << std::right << "0" << ": r_norm = " << std::fixed << std::setw(20) << std::setprecision(8) << std::right << r0_norm << " rel_resid = " << std::setw(14) << std::setprecision(10) << std::right << r0_norm/r0_init_norm << endl; } } @@ -252,7 +256,7 @@ void CostFunction::truncatedNewton(real* qstate, real* g, const real ftol) //CG cconvergence check r_norm = sqrt(rr); - rel_resid = r_norm/r0_norm; + rel_resid = r_norm/r0_init_norm; std::cout << std::right; if (verbose) cout << "\t\tCG iteration " << std::setw(7) << std::right << cg_its + 1 << ": r_norm = " << std::setw(20) << std::fixed << std::setprecision(8) << std::right << r_norm << " rel_resid = " << std::setw(14) << std::setprecision(10) << std::right << rel_resid << endl; diff --git a/src/CostFunction3D.cpp b/src/CostFunction3D.cpp index aa1ca33..b7eeb30 100644 --- a/src/CostFunction3D.cpp +++ b/src/CostFunction3D.cpp @@ -137,65 +137,91 @@ void CostFunction3D::finalize() void CostFunction3D::initialize(HashMap* config, real* bgU, real* obs, ReferenceState* ref, - uint64_t numVar, int numDerivatives, int numObMetaData) + uint64_t numVar, int numDerivatives, int numObMetaData,int aMode) { // Initialize number of variables - varDim = numVar; - derivDim = numDerivatives; - obMetaSize = numObMetaData; - configHash = config; + varDim = numVar; + analysisMode = aMode; + derivDim = numDerivatives; + obMetaSize = numObMetaData; + configHash = config; /* Set the output path */ dataPath = (*configHash)["data_directory"]; outputPath = (*configHash)["output_directory"]; - //JMD TODO: Need to plug in the Thermo boundary conditions here. - //JMD Currently this is only the wind boundary condidtions. // Horizontal boundary conditions - iBCL[0] = bcHash[(*configHash)["i_rhou_bcL"]]; - iBCR[0] = bcHash[(*configHash)["i_rhou_bcR"]]; - iBCL[1] = bcHash[(*configHash)["i_rhov_bcL"]]; - iBCR[1] = bcHash[(*configHash)["i_rhov_bcR"]]; - iBCL[2] = bcHash[(*configHash)["i_rhow_bcL"]]; - iBCR[2] = bcHash[(*configHash)["i_rhow_bcR"]]; - iBCL[3] = bcHash[(*configHash)["i_tempk_bcL"]]; - iBCR[3] = bcHash[(*configHash)["i_tempk_bcR"]]; - iBCL[4] = bcHash[(*configHash)["i_qv_bcL"]]; - iBCR[4] = bcHash[(*configHash)["i_qv_bcR"]]; - iBCL[5] = bcHash[(*configHash)["i_rhoa_bcL"]]; - iBCR[5] = bcHash[(*configHash)["i_rhoa_bcR"]]; - iBCL[6] = bcHash[(*configHash)["i_qr_bcL"]]; - iBCR[6] = bcHash[(*configHash)["i_qr_bcR"]]; - - jBCL[0] = bcHash[(*configHash)["j_rhou_bcL"]]; - jBCR[0] = bcHash[(*configHash)["j_rhou_bcR"]]; - jBCL[1] = bcHash[(*configHash)["j_rhov_bcL"]]; - jBCR[1] = bcHash[(*configHash)["j_rhov_bcR"]]; - jBCL[2] = bcHash[(*configHash)["j_rhow_bcL"]]; - jBCR[2] = bcHash[(*configHash)["j_rhow_bcR"]]; - jBCL[3] = bcHash[(*configHash)["j_tempk_bcL"]]; - jBCR[3] = bcHash[(*configHash)["j_tempk_bcR"]]; - jBCL[4] = bcHash[(*configHash)["j_qv_bcL"]]; - jBCR[4] = bcHash[(*configHash)["j_qv_bcR"]]; - jBCL[5] = bcHash[(*configHash)["j_rhoa_bcL"]]; - jBCR[5] = bcHash[(*configHash)["j_rhoa_bcR"]]; - jBCL[6] = bcHash[(*configHash)["j_qr_bcL"]]; - jBCR[6] = bcHash[(*configHash)["j_qr_bcR"]]; - - kBCL[0] = bcHash[(*configHash)["k_rhou_bcL"]]; - kBCR[0] = bcHash[(*configHash)["k_rhou_bcR"]]; - kBCL[1] = bcHash[(*configHash)["k_rhov_bcL"]]; - kBCR[1] = bcHash[(*configHash)["k_rhov_bcR"]]; - kBCL[2] = bcHash[(*configHash)["k_rhow_bcL"]]; - kBCR[2] = bcHash[(*configHash)["k_rhow_bcR"]]; - kBCL[3] = bcHash[(*configHash)["k_tempk_bcL"]]; - kBCR[3] = bcHash[(*configHash)["k_tempk_bcR"]]; - kBCL[4] = bcHash[(*configHash)["k_qv_bcL"]]; - kBCR[4] = bcHash[(*configHash)["k_qv_bcR"]]; - kBCL[5] = bcHash[(*configHash)["k_rhoa_bcL"]]; - kBCR[5] = bcHash[(*configHash)["k_rhoa_bcR"]]; - kBCL[6] = bcHash[(*configHash)["k_qr_bcL"]]; - kBCR[6] = bcHash[(*configHash)["k_qr_bcR"]]; + //if ((*configHash)["analysis_type"] == "WIND") { + if (analysisMode == VarDriver::WIND) { + std::cout << "Horizontal boundary condidtions for WIND" << std::endl; + iBCL[0] = bcHash[(*configHash)["i_rhou_bcL"]]; + iBCR[0] = bcHash[(*configHash)["i_rhou_bcR"]]; + iBCL[1] = bcHash[(*configHash)["i_rhov_bcL"]]; + iBCR[1] = bcHash[(*configHash)["i_rhov_bcR"]]; + iBCL[2] = bcHash[(*configHash)["i_rhow_bcL"]]; + iBCR[2] = bcHash[(*configHash)["i_rhow_bcR"]]; + iBCL[3] = bcHash[(*configHash)["i_tempk_bcL"]]; + iBCR[3] = bcHash[(*configHash)["i_tempk_bcR"]]; + iBCL[4] = bcHash[(*configHash)["i_qv_bcL"]]; + iBCR[4] = bcHash[(*configHash)["i_qv_bcR"]]; + iBCL[5] = bcHash[(*configHash)["i_rhoa_bcL"]]; + iBCR[5] = bcHash[(*configHash)["i_rhoa_bcR"]]; + iBCL[6] = bcHash[(*configHash)["i_qr_bcL"]]; + iBCR[6] = bcHash[(*configHash)["i_qr_bcR"]]; + + jBCL[0] = bcHash[(*configHash)["j_rhou_bcL"]]; + jBCR[0] = bcHash[(*configHash)["j_rhou_bcR"]]; + jBCL[1] = bcHash[(*configHash)["j_rhov_bcL"]]; + jBCR[1] = bcHash[(*configHash)["j_rhov_bcR"]]; + jBCL[2] = bcHash[(*configHash)["j_rhow_bcL"]]; + jBCR[2] = bcHash[(*configHash)["j_rhow_bcR"]]; + jBCL[3] = bcHash[(*configHash)["j_tempk_bcL"]]; + jBCR[3] = bcHash[(*configHash)["j_tempk_bcR"]]; + jBCL[4] = bcHash[(*configHash)["j_qv_bcL"]]; + jBCR[4] = bcHash[(*configHash)["j_qv_bcR"]]; + jBCL[5] = bcHash[(*configHash)["j_rhoa_bcL"]]; + jBCR[5] = bcHash[(*configHash)["j_rhoa_bcR"]]; + jBCL[6] = bcHash[(*configHash)["j_qr_bcL"]]; + jBCR[6] = bcHash[(*configHash)["j_qr_bcR"]]; + + kBCL[0] = bcHash[(*configHash)["k_rhou_bcL"]]; + kBCR[0] = bcHash[(*configHash)["k_rhou_bcR"]]; + kBCL[1] = bcHash[(*configHash)["k_rhov_bcL"]]; + kBCR[1] = bcHash[(*configHash)["k_rhov_bcR"]]; + kBCL[2] = bcHash[(*configHash)["k_rhow_bcL"]]; + kBCR[2] = bcHash[(*configHash)["k_rhow_bcR"]]; + kBCL[3] = bcHash[(*configHash)["k_tempk_bcL"]]; + kBCR[3] = bcHash[(*configHash)["k_tempk_bcR"]]; + kBCL[4] = bcHash[(*configHash)["k_qv_bcL"]]; + kBCR[4] = bcHash[(*configHash)["k_qv_bcR"]]; + kBCL[5] = bcHash[(*configHash)["k_rhoa_bcL"]]; + kBCR[5] = bcHash[(*configHash)["k_rhoa_bcR"]]; + kBCL[6] = bcHash[(*configHash)["k_qr_bcL"]]; + kBCR[6] = bcHash[(*configHash)["k_qr_bcR"]]; + //} else if((*configHash)["analysis_type"] == "THERMO") { + } else if(analysisMode == VarDriver::THERMO) { + std::cout << "Horizontal boundary condidtions for THERMO" << std::endl; + iBCL[0] = bcHash[(*configHash)["i_pip_bcL"]]; + iBCR[0] = bcHash[(*configHash)["i_pip_bcR"]]; + iBCL[1] = bcHash[(*configHash)["i_thetarhop_bcL"]]; + iBCR[1] = bcHash[(*configHash)["i_thetarhop_bcR"]]; + iBCL[2] = bcHash[(*configHash)["i_ftheta_bcL"]]; + iBCR[2] = bcHash[(*configHash)["i_ftheta_bcR"]]; + + jBCL[0] = bcHash[(*configHash)["j_pip_bcL"]]; + jBCR[0] = bcHash[(*configHash)["j_pip_bcR"]]; + jBCL[1] = bcHash[(*configHash)["j_thetarhop_bcL"]]; + jBCR[1] = bcHash[(*configHash)["j_thetarhop_bcR"]]; + jBCL[2] = bcHash[(*configHash)["j_ftheta_bcL"]]; + jBCR[2] = bcHash[(*configHash)["j_ftheta_bcR"]]; + + kBCL[0] = bcHash[(*configHash)["k_pip_bcL"]]; + kBCR[0] = bcHash[(*configHash)["k_pip_bcR"]]; + kBCL[1] = bcHash[(*configHash)["k_thetarhop_bcL"]]; + kBCR[1] = bcHash[(*configHash)["k_thetarhop_bcR"]]; + kBCL[2] = bcHash[(*configHash)["k_ftheta_bcL"]]; + kBCR[2] = bcHash[(*configHash)["k_ftheta_bcR"]]; + } // Define the Reference state refstate = ref; @@ -250,6 +276,7 @@ void CostFunction3D::initialize(HashMap* config, CTHTd = new real[nState]; stateU = new real[nState]; int64_t vector_size = mObs*(obMetaSize+varDim*derivDim); + std::cout << "CostFunction3D::initialize: " << mObs << " " << vector_size << std::endl; obsVector = new real[vector_size]; obsData = new real[mObs]; HCq = new real[mObs+nodes]; @@ -278,11 +305,11 @@ void CostFunction3D::initialize(HashMap* config, kRankMax=0; for (int var = 0; var < varDim; ++var) { iRank[var] = iDim - rankHash[iBCL[var]] - rankHash[iBCR[var]]; - if (iBCL[var] == PERIODIC) {cout << "Periodic in I dimension\n"; iRank[var]--;} + if (iBCL[var] == PERIODIC) {std::cout << "Periodic in I dimension" << std::endl; iRank[var]--;} jRank[var] = jDim - rankHash[jBCL[var]] - rankHash[jBCR[var]]; - if (jBCL[var] == PERIODIC) {cout << "Periodic in J dimension\n"; jRank[var]--;} + if (jBCL[var] == PERIODIC) {std::cout << "Periodic in J dimension" << std::endl; jRank[var]--;} kRank[var] = kDim - rankHash[kBCL[var]] - rankHash[kBCR[var]]; - if (kBCL[var] == PERIODIC) {cout << "Periodic in K dimension\n"; kRank[var]--;} + if (kBCL[var] == PERIODIC) {std::cout << "Periodic in K dimension" << std::endl; kRank[var]--;} // Need to verify that Rank is sufficient (>2?) iL[var] = new real[iRank[var]*iLDim]; jL[var] = new real[jRank[var]*jLDim]; @@ -293,7 +320,7 @@ void CostFunction3D::initialize(HashMap* config, kRankMax = max(kRank[var],kRankMax); } - cout << "kRankMax: " << kRankMax << "\n"; + std::cout << "kRankMax: " << kRankMax << std::endl; /* Precalculate the basis functions for lookup table option basisappx = configHash->value("spline_approximation").toInt(); @@ -327,7 +354,7 @@ void CostFunction3D::initState(const int iteration) { GPTLstart("CostFunction3D::initState"); // Clear the state vector - cout << "Initializing state vector..." << endl; + std::cout << "Initializing state vector..." << std::endl; for (int n = 0; n < nState; n++) { currState[n] = 0.0; currGradient[n] = 0.0; @@ -355,7 +382,9 @@ void CostFunction3D::initState(const int iteration) jMaxWavenumber[i] = -1.0; kMaxWavenumber[i] = -1.0; } - if (configHash->exists("i_max_wavenumber")) { + // JMD Not quite sure what is going on here. In the TDRP files only + // JMD three values exist yet this appears to be working for 7 variables + if (configHash->exists("i_max_wavenumber")) { // Set all the variables to the same filter for (int i = 0; i < varDim; i++) { iMaxWavenumber[i] = std::stof((*configHash)["i_max_wavenumber"]); @@ -371,7 +400,7 @@ void CostFunction3D::initState(const int iteration) iMaxWavenumber[6] = std::stof((*configHash)["i_max_wavenumber_qr"]); } - if (configHash->exists("j_max_wavenumber")) { + if (configHash->exists("j_max_wavenumber")) { for (int i = 0; i < varDim; i++) { jMaxWavenumber[i] = std::stof((*configHash)["j_max_wavenumber"]); } @@ -385,7 +414,7 @@ void CostFunction3D::initState(const int iteration) jMaxWavenumber[6] = std::stof((*configHash)["j_max_wavenumber_qr"]); } - if (configHash->exists("k_max_wavenumber")) { + if (configHash->exists("k_max_wavenumber")) { for (int i = 0; i < varDim; i++) { kMaxWavenumber[i] = std::stof((*configHash)["k_max_wavenumber"]); } @@ -403,12 +432,12 @@ void CostFunction3D::initState(const int iteration) UseFFT=false; for (int var = 0; var < varDim; var++) { if ((kBCL[var] == PERIODIC) and (kMaxWavenumber[var] >= 0)) UseFFT=true; - // cout << "jBCL[var]: " << jBCL[var] << " jMaxWavenumber[var]: " << jMaxWavenumber[var] << "\n"; + // std::cout << "jBCL[var]: " << jBCL[var] << " jMaxWavenumber[var]: " << jMaxWavenumber[var] << std::endl; if ((jBCL[var] == PERIODIC) and (jMaxWavenumber[var] >= 0)) UseFFT=true; if ((iBCL[var] == PERIODIC) and (iMaxWavenumber[var] >= 0)) UseFFT=true; } if(UseFFT) - cout << "PERIODIC boundaries and maximum wavenumber enforcement will enable FFtransform() \n"; + std::cout << "PERIODIC boundaries and maximum wavenumber enforcement will enable FFtransform()" << std::endl; // Set up the spline matrices @@ -424,11 +453,11 @@ void CostFunction3D::initState(const int iteration) // Mass continuity weight mcWeight = std::stof((*configHash)["mc_weight"]); - cout << "Mass continuity weight set to " << mcWeight << endl; + std::cout << "Mass continuity weight set to " << mcWeight << std::endl; if (iteration == 1) { - cout << "Initializing background..." << endl; + std::cout << "Initializing background..." << std::endl; // Set up the background state for (int n = 0; n < nState; n++) { bgState[n] = 0.0; @@ -450,9 +479,9 @@ void CostFunction3D::initState(const int iteration) // variance.writeDebugNc("debug.out/std_errors.nc", true, mishError); SBtransform(mishError, SBError); // variance.writeDebugNc("debug.out/std_errors_SB.nc", false, SBError); - #pragma acc data copyin(SBError[0:nState]) copyout(meshError[0:nState]) - { - SAtransform(SBError, meshError); + #pragma acc data copyin(SBError[0:nState]) copyout(meshError[0:nState]) + { + SAtransform(SBError, meshError); } variance.setMeshData(meshError); delete[] SBError; @@ -482,7 +511,7 @@ void CostFunction3D::initState(const int iteration) // FF transform to match background and increment #pragma acc data copyin(stateA[0:nState]) copyout(bgState[0:nState]) { - FFtransform(stateA, bgState); + FFtransform(stateA, bgState); } // variance.writeDebugNc("debug.out/FF.nc", false, bgState); // mesh sized DEBUG @@ -493,12 +522,10 @@ void CostFunction3D::initState(const int iteration) // Init node variance for (int iIndex = 0; iIndex < iDim; iIndex++) { for (int jIndex = 0; jIndex < jDim; jIndex++) { - for (int kIndex = 0; kIndex < kDim; kIndex++) { -// int bIndex = varDim * iDim * jDim*kIndex + varDim * iDim * jIndex + varDim * iIndex + var; - int64_t bIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); -// *bIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); - bgStdDev[bIndex] = variance.meshValueAt(var, iIndex, jIndex, kIndex); - } + for (int kIndex = 0; kIndex < kDim; kIndex++) { + int64_t bIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); + bgStdDev[bIndex] = variance.meshValueAt(var, iIndex, jIndex, kIndex); + } } } } @@ -509,30 +536,29 @@ void CostFunction3D::initState(const int iteration) real varScale = 0; for (int iIndex = 0; iIndex < iDim; iIndex++) { for (int jIndex = 0; jIndex < jDim; jIndex++) { - for (int kIndex = 0; kIndex < kDim; kIndex++) { -// int bIndex = varDim * iDim * jDim * kIndex + varDim * iDim * jIndex + varDim * iIndex + var; - int64_t bIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); -// int64_t *bIndex = (int64_t *)malloc(sizeof(int64_t)); -// *bIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); - varScale += bgState[bIndex] * bgState[bIndex]; - } + for (int kIndex = 0; kIndex < kDim; kIndex++) { + int64_t bIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); + varScale += bgState[bIndex] * bgState[bIndex]; + } } } varScale = sqrt(varScale / (iDim * jDim * kDim)); if (varScale) { real errPct = 100 * bgError[var] / varScale; - cout << "Variable " << var << " RMS = " << varScale << "\t BG Error = " << bgError[var] - << " ( " << errPct << " %)" << endl; + std::cout << "Variable " << var << " RMS = " << varScale << "\t BG Error = " << bgError[var] + << " ( " << errPct << " %)" << std::endl; } else { - cout << "Variable " << var << " RMS = " << varScale << "\t BG Error = " << bgError[var] - << " ( Infinite! %)" << endl; + std::cout << "Variable " << var << " RMS = " << varScale << "\t BG Error = " << bgError[var] + << " ( Infinite! %)" << std::endl; } } // Load the obs locally and weight the nonlinear observation operators by interpolated bg fields + std::cout << "Before call to obAdjustments" << std::endl; obAdjustments(); + std::cout << "Before call to obAdjustments" << std::endl; // Calculate the H matrix operator calcHmatrix(); @@ -540,24 +566,28 @@ void CostFunction3D::initState(const int iteration) calcInnovation(); // Output the original background field - outputAnalysis("background", bgState); + if(analysisMode == VarDriver::WIND) { + outputAnalysis("background", bgState); + } else if (analysisMode == VarDriver::THERMO) { + outputAnalysis_thermo("background", bgState); + } + - cout << "Beginning analysis...\n"; + std::cout << "Beginning analysis..." << std::endl; // HTd calcHTranspose(innovation, stateC); -#pragma acc data create(stateB[0:nState]) -{ - // cout << "Before FFtransform ...\n"; - FFtransform(stateC, stateA); - // cout << "Before SAtransform ...\n"; - SAtransform(stateA, stateB); - // cout << "Before SCtransform ...\n"; - SCtransform(stateB, CTHTd); -} + #pragma acc data create(stateB[0:nState]) + { + // std::cout << "Before FFtransform ..." << std::endl; + FFtransform(stateC, stateA); + // std::cout << "Before SAtransform ..." << std::endl; + SAtransform(stateA, stateB); + // std::cout << "Before SCtransform ..." << std::endl; + SCtransform(stateB, CTHTd); + } #pragma acc enter data copyin(CTHTd) - //Htransform(stateB); GPTLstop("CostFunction3D::initState"); } @@ -572,34 +602,34 @@ real CostFunction3D::funcValue(real* state) obIP = 0.; GPTLstart("CostFunction3D::funcValue"); - #pragma acc data copyin(state[0:nState]) - { + #pragma acc data copyin(state[0:nState]) + { - GPTLstart("CostFunction3D::funcValue:updateHCq"); - updateHCq(state,HCq); - GPTLstop("CostFunction3D::funcValue:updateHCq"); + GPTLstart("CostFunction3D::funcValue:updateHCq"); + updateHCq(state,HCq); + GPTLstop("CostFunction3D::funcValue:updateHCq"); - GPTLstart("CostFunction3D::funcValue:other"); - // Compute inner product of state vector - #pragma omp parallel for reduction(+:qIP) - #pragma acc parallel loop reduction(+:qIP) - for (int n = 0; n < nState; n++) { - qIP += state[n]*state[n]; - } + GPTLstart("CostFunction3D::funcValue:other"); + // Compute inner product of state vector + #pragma omp parallel for reduction(+:qIP) + #pragma acc parallel loop reduction(+:qIP) + for (int n = 0; n < nState; n++) { + qIP += state[n]*state[n]; + } - // Subtract d from HCq to yield mObs length vector and compute inner product - #pragma omp parallel for reduction(+:obIP) - #pragma acc parallel loop reduction(+:obIP) private(obIndex,m) - for (m = 0; m < mObs; m++) { - //obIP += (HCq[m]-innovation[m])*(obsVector[obIndex])*(HCq[m]-innovation[m]); - obIP += (HCq[m]-innovation[m])*(obsData[m])*(HCq[m]-innovation[m]); - } - GPTLstop("CostFunction3D::funcValue:other"); - } + // Subtract d from HCq to yield mObs length vector and compute inner product + #pragma omp parallel for reduction(+:obIP) + #pragma acc parallel loop reduction(+:obIP) private(obIndex,m) + for (m = 0; m < mObs; m++) { + //obIP += (HCq[m]-innovation[m])*(obsVector[obIndex])*(HCq[m]-innovation[m]); + obIP += (HCq[m]-innovation[m])*(obsData[m])*(HCq[m]-innovation[m]); + } + GPTLstop("CostFunction3D::funcValue:other"); + } - J = 0.5*(qIP + obIP); - GPTLstop("CostFunction3D::funcValue"); - return J; + J = 0.5*(qIP + obIP); + GPTLstop("CostFunction3D::funcValue"); + return J; } void CostFunction3D::funcGradient(real* state, real* gradient) @@ -607,39 +637,37 @@ void CostFunction3D::funcGradient(real* state, real* gradient) int n; GPTLstart("CostFunction3D::funcGradient"); - #pragma acc data copyin(state[0:nState]) copyout(gradient[0:nState]) - { + #pragma acc data copyin(state[0:nState]) copyout(gradient[0:nState]) + { + GPTLstart("CostFunction3D::funcGradient:updateHCq"); + updateHCq(state,HCq); + GPTLstop("CostFunction3D::funcGradient:updateHCq"); - GPTLstart("CostFunction3D::funcGradient:updateHCq"); - updateHCq(state,HCq); - GPTLstop("CostFunction3D::funcGradient:updateHCq"); - - GPTLstart("CostFunction3D::funcGradient:calcHTranspose"); - // HTHCq - calcHTranspose(HCq, stateC); - GPTLstop("CostFunction3D::funcGradient:calcHTranspose"); - - GPTLstart("CostFunction3D::funcGradient:FFtransform"); - FFtransform(stateC, stateA); - GPTLstop("CostFunction3D::funcGradient:FFtransform"); - - GPTLstart("CostFunction3D::funcGradient:SAtransform"); - SAtransform(stateA, stateB); - GPTLstop("CostFunction3D::funcGradient:SAtransform"); - GPTLstart("CostFunction3D::funcGradient:SCtransform"); - SCtransform(stateB, stateC); - GPTLstop("CostFunction3D::funcGradient:SCtransform"); - - GPTLstart("CostFunction3D::funcGradient:gradient"); - #pragma acc parallel loop gang vector vector_length(32) private(n) - for (n = 0; n < nState; n++) { - gradient[n] = state[n] + stateC[n] - CTHTd[n]; - } - GPTLstop("CostFunction3D::funcGradient:gradient"); + GPTLstart("CostFunction3D::funcGradient:calcHTranspose"); + // HTHCq + calcHTranspose(HCq, stateC); + GPTLstop("CostFunction3D::funcGradient:calcHTranspose"); - } + GPTLstart("CostFunction3D::funcGradient:FFtransform"); + FFtransform(stateC, stateA); + GPTLstop("CostFunction3D::funcGradient:FFtransform"); - GPTLstop("CostFunction3D::funcGradient"); + GPTLstart("CostFunction3D::funcGradient:SAtransform"); + SAtransform(stateA, stateB); + GPTLstop("CostFunction3D::funcGradient:SAtransform"); + GPTLstart("CostFunction3D::funcGradient:SCtransform"); + SCtransform(stateB, stateC); + GPTLstop("CostFunction3D::funcGradient:SCtransform"); + + GPTLstart("CostFunction3D::funcGradient:gradient"); + #pragma acc parallel loop gang vector vector_length(32) private(n) + for (n = 0; n < nState; n++) { + gradient[n] = state[n] + stateC[n] - CTHTd[n]; + } + GPTLstop("CostFunction3D::funcGradient:gradient"); + + } + GPTLstop("CostFunction3D::funcGradient"); } @@ -655,46 +683,45 @@ real CostFunction3D::funcValueAndGradient(real* state, real *gradient) #pragma acc data present(state[0:nState],gradient[0:nState]) { - qIP = 0.; - obIP = 0.; - // cout << "CostFunction3D::funcValueAndGradient\n"; - - //update global HCq variable - updateHCq(state,HCq); + qIP = 0.; + obIP = 0.; + + //update global HCq variable + updateHCq(state,HCq); + + //Func Value + // Compute inner product of state vector + //#pragma omp parallel for reduction(+:qIP) + #pragma acc parallel loop reduction(+:qIP) + for (n = 0; n < nState; n++) { + qIP += state[n]*state[n]; + } + // Subtract d from HCq to yield mObs length vector and compute inner product + //#pragma omp parallel for reduction(+:obIP) + #pragma acc parallel loop reduction(+:obIP) private(m) + for (m = 0; m < mObs; m++) { + //int obIndex = m*(obMetaSize+varDim*derivDim) + 1; + //obIP += (HCq[m]-innovation[m])*(obsVector[obIndex])*(HCq[m]-innovation[m]); + obIP += (HCq[m]-innovation[m])*(obsData[m])*(HCq[m]-innovation[m]); + } + //function value J + J = 0.5*(qIP + obIP); - //Func Value - // Compute inner product of state vector - //#pragma omp parallel for reduction(+:qIP) - #pragma acc parallel loop reduction(+:qIP) - for (n = 0; n < nState; n++) { - qIP += state[n]*state[n]; - } - // Subtract d from HCq to yield mObs length vector and compute inner product - //#pragma omp parallel for reduction(+:obIP) - #pragma acc parallel loop reduction(+:obIP) private(m) - for (m = 0; m < mObs; m++) { - //int obIndex = m*(obMetaSize+varDim*derivDim) + 1; - //obIP += (HCq[m]-innovation[m])*(obsVector[obIndex])*(HCq[m]-innovation[m]); - obIP += (HCq[m]-innovation[m])*(obsData[m])*(HCq[m]-innovation[m]); - } - //function value J - J = 0.5*(qIP + obIP); - - //Now Gradient (also uses HCq) - calcHTranspose(HCq, stateC); - FFtransform(stateC, stateA); - SAtransform(stateA, stateB); - SCtransform(stateB, stateC); - //calc gradient - #pragma acc parallel loop gang vector vector_length(32) - for (int n = 0; n < nState; n++) { - gradient[n] = state[n] + stateC[n] - CTHTd[n]; - } + //Now Gradient (also uses HCq) + calcHTranspose(HCq, stateC); + FFtransform(stateC, stateA); + SAtransform(stateA, stateB); + SCtransform(stateB, stateC); + //calc gradient + #pragma acc parallel loop gang vector vector_length(32) + for (int n = 0; n < nState; n++) { + gradient[n] = state[n] + stateC[n] - CTHTd[n]; + } - GPTLstop("CostFunction3D::funcValueAndGradient"); - //cout << "CostFunction3D::funcValueAndGradient: qIP, obIP " << qIP << " " << obIP << "\n"; + GPTLstop("CostFunction3D::funcValueAndGradient"); + //std::cout << "CostFunction3D::funcValueAndGradient: qIP, obIP " << qIP << " " << obIP << std::endl; } - return J; + return J; } @@ -705,7 +732,7 @@ void CostFunction3D::funcHessian(real* x, real *hessian) int n; #pragma acc data present(x[0:nState],hessian[0:nState]) { - //DBG cout << "CostFunction3D::funcHessian\n"; + //DBG std::cout << "CostFunction3D::funcHessian" << std::endl; //calc HCx (store in global variable HCq) updateHCq(x,HCq); @@ -747,35 +774,44 @@ void CostFunction3D::updateBG() // S (SA transform) yield A's #pragma acc data copyin(currState[0:nState]) copyout(stateC[0:nState]) create(stateB[0:nState]) { - SCtransform(currState, stateB); - SAtransform(stateB, stateA); - FFtransform(stateA, stateC); + SCtransform(currState, stateB); + SAtransform(stateB, stateA); + FFtransform(stateA, stateC); } - outputAnalysis("increment", stateC); + if(analysisMode == VarDriver::WIND) { + outputAnalysis("increment", stateC); + } else if (analysisMode == VarDriver::THERMO) { + outputAnalysis_thermo("increment", stateC); + } + // In BG update we are directly summing C + A std::string cFilename = outputPath + "/samurai_Coefficients.out"; ofstream cstream(cFilename); //JMD variable-interleave - cstream << "Variable\tI\tJ\tK\tBackground\tAnalysis\tIncrement\n"; + cstream << "Variable\tI\tJ\tK\tBackground\tAnalysis\tIncrement" << std::endl; for (int var = 0; var < varDim; var++) { for (int iIndex = 0; iIndex < iDim; iIndex++) { for (int jIndex = 0; jIndex < jDim; jIndex++) { - for (int kIndex = 0; kIndex < kDim; kIndex++) { - cstream << var << "\t" << iIndex << "\t" << jIndex << "\t" << kIndex << "\t"; - //int bgIndex = varDim*iDim*jDim*kIndex + varDim*iDim*jIndex +varDim*iIndex + var; - int bgIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); - cstream << bgState[bgIndex] << "\t"; - bgState[bgIndex] += stateC[bgIndex]; - cstream << bgState[bgIndex] << "\t"; - cstream << stateC[bgIndex] << endl; - } + for (int kIndex = 0; kIndex < kDim; kIndex++) { + cstream << var << "\t" << iIndex << "\t" << jIndex << "\t" << kIndex << "\t"; + int bgIndex = INDEX(iIndex, jIndex, kIndex, iDim, jDim, varDim, var); + cstream << bgState[bgIndex] << "\t"; + bgState[bgIndex] += stateC[bgIndex]; + cstream << bgState[bgIndex] << "\t"; + cstream << stateC[bgIndex] << std::endl; + } } } } - outputAnalysis("analysis", bgState); + if(analysisMode == VarDriver::WIND) { + outputAnalysis("analysis", bgState); + } else if (analysisMode == VarDriver::THERMO) { + outputAnalysis_thermo("analysis", bgState); + } + GPTLstop("CostFunction3D::updateBG"); } @@ -783,13 +819,13 @@ void CostFunction3D::calcInnovation() { GPTLstart("CostFunction3D::calcInnovation"); // Initialize and fill the innovation vector - cout << "Initializing innovation vector..." << endl; + std::cout << "Initializing innovation vector..." << std::endl; // Use HCq to hold the transform, but C is not applied for the innovation -#pragma acc data copyin(bgState[0:nState]) create(HCq[mObs+(iDim*jDim*kDim)]) -{ - Htransform(bgState, HCq); -} + #pragma acc data copyin(bgState[0:nState]) create(HCq[mObs+(iDim*jDim*kDim)]) + { + Htransform(bgState, HCq); + } real innovationRMS = 0.; @@ -804,7 +840,7 @@ void CostFunction3D::calcInnovation() if (mObs) innovationRMS /= mObs; innovationRMS = sqrt(innovationRMS); - cout << "Innovation RMS : " << innovationRMS << endl; + std::cout << "Innovation RMS : " << innovationRMS << std::endl; #pragma acc enter data copyin(innovation) GPTLstop("CostFunction3D::calcInnovation"); } @@ -988,7 +1024,7 @@ bool CostFunction3D::SAtransform(const real* Bstate, real* Astate) for (int m = 0; m < iRank[var]; m++) { tmp += iGamma[var][iDim*m + i]*xi[m]; } - // std::cout << "i: " << i << " ai[" << i << "]: " << tmp << "\n"; + // std::cout << "i: " << i << " ai[" << i << "]: " << tmp << std::endl; Astate[INDEX(i, jIndex, kIndex, iDim, jDim, varDim, var)] = tmp; } } @@ -1128,7 +1164,7 @@ bool CostFunction3D::SAtranspose(const real* Astate, real* Bstate) for (int k = 0; k < kDim; k++) { b[m] += kGamma[var][kDim*m + k]*kB[k]; } - //std::cout << m << " " << b[m] << "\n"; + //std::cout << m << " " << b[m] << std::endl; } // Solve for A's using compact storage @@ -1154,7 +1190,7 @@ bool CostFunction3D::SAtranspose(const real* Astate, real* Bstate) for (int m = 0; m < kRank[var]; m++) { a[k] += kGamma[var][kDim*m + k]*x[m]; } - //std::cout << k << " " << a[k] << "\n"; + //std::cout << k << " " << a[k] << std::endl; } for (int kIndex = 0; kIndex < kDim; kIndex++) { @@ -1197,7 +1233,7 @@ void CostFunction3D::SBtransform(const real* Ustate, real* Bstate) real gausspoint = 0.5*sqrt(1./3.); for (var = 0; var < varDim; var++) { - cout << " min(rankHash[iBCL[var]],1): " << min(rankHash[iBCL[var]],1) << " max(iDim-1-rankHash[iBCR[var]],iDim-2): " << max(iDim-1-rankHash[iBCR[var]],iDim-2) << "\n"; + std::cout << " min(rankHash[iBCL[var]],1): " << min(rankHash[iBCL[var]],1) << " max(iDim-1-rankHash[iBCR[var]],iDim-2): " << max(iDim-1-rankHash[iBCR[var]],iDim-2) << std::endl; is = min(rankHash[iBCL[var]],1); ie = max(iDim-1-rankHash[iBCR[var]],iDim-2); js = min(rankHash[jBCL[var]],1); @@ -1624,17 +1660,17 @@ bool CostFunction3D::setupSplines() real Pi = acos(-1.); real cutoff_wl = std::stof((*configHash)["i_spline_cutoff"]); - cout << "i Spline cutoff set to " << cutoff_wl << endl; + std::cout << "i Spline cutoff set to " << cutoff_wl << std::endl; real eq = pow( (cutoff_wl/(2*Pi)) , 6); calcSplineCoefficients(iDim, eq, iBCL, iBCR, iMin, DI, DIrecip, iLDim, iL, iGamma); cutoff_wl = std::stof((*configHash)["j_spline_cutoff"]); - cout << "j Spline cutoff set to " << cutoff_wl << endl; + std::cout << "j Spline cutoff set to " << cutoff_wl << std::endl; eq = pow( (cutoff_wl/(2*Pi)) , 6); calcSplineCoefficients(jDim, eq, jBCL, jBCR, jMin, DJ, DJrecip, jLDim, jL, jGamma); cutoff_wl = std::stof((*configHash)["k_spline_cutoff"]); - cout << "k Spline cutoff set to " << cutoff_wl << endl; + std::cout << "k Spline cutoff set to " << cutoff_wl << std::endl; eq = pow( (cutoff_wl/(2*Pi)) , 6); calcSplineCoefficients(kDim, eq, kBCL, kBCR, kMin, DK, DKrecip, kLDim, kL, kGamma); @@ -1647,6 +1683,7 @@ bool CostFunction3D::setupSplines() void CostFunction3D::obAdjustments() { GPTLstart("CostFunction3D::obAdjustments"); + std::cout << "CostFunction3D::obAdjustments: mObs " << mObs << std::endl; //JMD variable-interleave // Load the obs locally and weight the nonlinear observation operators by interpolated bg fields for (int m = 0; m < mObs; m++) { @@ -1665,8 +1702,8 @@ void CostFunction3D::obAdjustments() { if ((i < iMin) or (i > iMax) or (j < jMin) or (j > jMax) or (k < kMin) or (k > kMax)) { - cout << "Error! Observations are found outside the domain where the spline is undefined.\n"; - cout << "This can only happen if you bypassed preprocessing -- check your samurai_Observations.in and re-run.\n"; + std::cout << "Error! Observations are found outside the domain where the spline is undefined." << std::endl; + std::cout << "This can only happen if you bypassed preprocessing -- check your samurai_Observations.in and re-run." << std::endl; } real rhoprime = 0.; real qvprime = 0.; @@ -1693,7 +1730,7 @@ void CostFunction3D::obAdjustments() { // *bIndex = varDim * iDim * jDim * kNode + varDim * iDim * jNode + varDim * iNode + 4; // *bIndex2 = varDim * iDim * jDim * kNode + varDim * iDim * jNode + varDim * iNode + 5; // turn off print out January 24, 2023 - // cout << "bIndex = " << *bIndex << endl; + // std::cout << "bIndex = " << *bIndex << std::endl; ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 0, iBCL[4], iBCR[4]); jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 0, jBCL[4], jBCR[4]); kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 0, kBCL[4], kBCR[4]); @@ -1731,6 +1768,7 @@ void CostFunction3D::obAdjustments() { for(int m=0;m= 0)) { - //cout << "FFT perform on J-pencil\n"; + //std::cout << "FFT perform on J-pencil" << std::endl; for (kIndex = 0; kIndex < kDim; kIndex++) { for (iIndex = 0; iIndex < iDim; iIndex++) { for (jIndex = 0; jIndex < jDim; jIndex++) { diff --git a/src/CostFunction3D.h b/src/CostFunction3D.h index df4a03d..58f71e6 100644 --- a/src/CostFunction3D.h +++ b/src/CostFunction3D.h @@ -42,7 +42,7 @@ class CostFunction3D : public CostFunction CostFunction3D(const Projection& proj, const int& numObs = 0, const int& stateSize = 0); virtual ~CostFunction3D(); - void initialize(HashMap* config, real* bgU, real* obs, ReferenceState* ref, uint64_t numVar, int numDerivatives, int numObMetaData); + void initialize(HashMap* config, real* bgU, real* obs, ReferenceState* ref, uint64_t numVar, int numDerivatives, int numObMetaData, int analysisMode); void finalize(); void updateBG(); void initState(const int iteration); @@ -72,6 +72,7 @@ CostFunction3D(const Projection& proj, const int& numObs = 0, const int& stateSi void calcInnovation(); void calcHTranspose(const real* yhat, real* Astate); virtual bool outputAnalysis(const std::string& suffix, real* Astate) = 0; + virtual bool outputAnalysis_thermo(const std::string& suffix, real* Astate) = 0; void SBtransform(const real* Ustate, real* Bstate); void SBtranspose(const real* Bstate, real* Ustate); void SCtransform(const real* Astate, real* Cstate); @@ -142,6 +143,7 @@ CostFunction3D(const Projection& proj, const int& numObs = 0, const int& stateSi uint64_t varDim; // NCAR: promoted to 64-bit, since it should auto-promote calculations with it to 64-bit int derivDim; int obMetaSize; + int analysisMode; real bgError[7]; int iBCL[7], iBCR[7], jBCL[7], jBCR[7], kBCL[7], kBCR[7]; diff --git a/src/CostFunctionCOAMPS.cpp b/src/CostFunctionCOAMPS.cpp index df03f7a..e56f22e 100644 --- a/src/CostFunctionCOAMPS.cpp +++ b/src/CostFunctionCOAMPS.cpp @@ -28,6 +28,9 @@ void CostFunctionCOAMPS::setSigmas(float *sigmas, int size) { sDim = size; } +bool CostFunctionCOAMPS::outputAnalysis_thermo(const std::string& suffix, real* Astate) +{ +} bool CostFunctionCOAMPS::outputAnalysis(const std::string& suffix, real* Astate) { /* NCAR FIXME diff --git a/src/CostFunctionCOAMPS.h b/src/CostFunctionCOAMPS.h index a2fcfeb..604e846 100644 --- a/src/CostFunctionCOAMPS.h +++ b/src/CostFunctionCOAMPS.h @@ -22,6 +22,7 @@ class CostFunctionCOAMPS: 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 writeFlatfile(const std::string& flatFileName, const int var); diff --git a/src/CostFunctionRTZ.cpp b/src/CostFunctionRTZ.cpp index ae3a396..1cb74bc 100644 --- a/src/CostFunctionRTZ.cpp +++ b/src/CostFunctionRTZ.cpp @@ -22,10 +22,13 @@ CostFunctionRTZ::CostFunctionRTZ(const Projection& proj, const int& numObs, cons CostFunctionRTZ::~CostFunctionRTZ() { } +bool CostFunctionRTZ::outputAnalysis_thermo(const std::string& suffix, real* Astate) +{ +} bool CostFunctionRTZ::outputAnalysis(const std::string& suffix, real* Astate) { std::string samuraiMode = "wind"; - cout << "Outputting " << suffix << "...\n"; + std::cout << "Outputting " << suffix << "...\n"; // H --> to Mish for output std::string samuraiout = "samurai_RTZ_" + samuraiMode + "_" + suffix + ".out"; ofstream samuraistream; @@ -542,13 +545,13 @@ bool CostFunctionRTZ::outputAnalysis(const std::string& suffix, real* Astate) if ((*configHash)["output_netcdf"] == "true") { std::string cdfFileName = outFileName + ".nc"; if (!writeNetCDF(cdfFileName)) - cout << "Error writing netcdf file " << cdfFileName << endl; + 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(asiFileName)) - cout << "Error writing asi file " << asiFileName << endl; + std::cout << "Error writing asi file " << asiFileName << std::endl; } // Set the domain back adjustInternalDomain(1); diff --git a/src/CostFunctionRTZ.h b/src/CostFunctionRTZ.h index 53f2cdf..c8c7f58 100644 --- a/src/CostFunctionRTZ.h +++ b/src/CostFunctionRTZ.h @@ -21,6 +21,7 @@ class CostFunctionRTZ: 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); }; diff --git a/src/CostFunctionXYP.cpp b/src/CostFunctionXYP.cpp index af1ca72..d3c426b 100644 --- a/src/CostFunctionXYP.cpp +++ b/src/CostFunctionXYP.cpp @@ -23,6 +23,10 @@ CostFunctionXYP::~CostFunctionXYP() { } +bool CostFunctionXYP::outputAnalysis_thermo(const std::string& suffix, real* Astate) +{ +} + bool CostFunctionXYP::outputAnalysis(const std::string& suffix, real* Astate) { // H --> to Mish for output diff --git a/src/CostFunctionXYP.h b/src/CostFunctionXYP.h index 16f3d7c..f5587aa 100644 --- a/src/CostFunctionXYP.h +++ b/src/CostFunctionXYP.h @@ -21,6 +21,7 @@ class CostFunctionXYP: 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); real pMin, pMax, DP; diff --git a/src/CostFunctionXYZ.cpp b/src/CostFunctionXYZ.cpp index ee632d7..53e7c5c 100644 --- a/src/CostFunctionXYZ.cpp +++ b/src/CostFunctionXYZ.cpp @@ -678,13 +678,13 @@ bool CostFunctionXYZ::outputAnalysis(const std::string& suffix, real* Astate) if ((*configHash)["output_netcdf"] == "true") { std::string cdfFileName = outFileName + ".nc"; if (!writeNetCDF(cdfFileName)) - cout << "Error writing netcdf file " << cdfFileName << endl; + 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(asiFileName)) - cout << "Error writing asi file " << asiFileName << endl; + std::cout << "Error writing asi file " << asiFileName << std::endl; } // Set the domain back adjustInternalDomain(1); @@ -827,7 +827,7 @@ bool CostFunctionXYZ::outputAnalysis_thermo(const std::string& suffix, real* Ast } std::string fileName = "samurai_XYZ_" + samuraiMode + "_" + suffix; - std::string outFileName = outputPath + fileName; + std::string outFileName = outputPath + "/" + fileName; // Write the Obs to a summary text file if ((*configHash)["output_qc"] == "true") { diff --git a/src/ErrorData.cpp b/src/ErrorData.cpp index b387fa2..f385ab5 100644 --- a/src/ErrorData.cpp +++ b/src/ErrorData.cpp @@ -10,9 +10,10 @@ double *ErrorData::init(std::string fname, HashMap* config, size_t numVar) { configHash = config; varDim = numVar; + + // JMD KLUDGE: Don't know how to address this code for the THERMO version // set defaults from the config hash - bgError[0] = std::stof((*configHash)["bg_rhou_error"]); bgError[1] = std::stof((*configHash)["bg_rhov_error"]); bgError[2] = std::stof((*configHash)["bg_rhow_error"]); diff --git a/src/NetCDF.h b/src/NetCDF.h index f638cd1..a41cd2d 100644 --- a/src/NetCDF.h +++ b/src/NetCDF.h @@ -23,8 +23,9 @@ class NetCDF virtual ~NetCDF(); - virtual int readNetCDF(const char* filename)=0; + virtual int readNetCDF(std::string filename)=0; virtual double getValue(const int &i,const int &j,const int &k,const std::string& varName)=0; + virtual int getValue(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; diff --git a/src/NetCDF_XYZ.cpp b/src/NetCDF_XYZ.cpp index eff84eb..dfdeeed 100644 --- a/src/NetCDF_XYZ.cpp +++ b/src/NetCDF_XYZ.cpp @@ -68,21 +68,29 @@ NetCDF_XYZ::~NetCDF_XYZ() } -int NetCDF_XYZ::readNetCDF(const char* filename) { +int NetCDF_XYZ::readNetCDF(std::string filename) { Nc3Error err(Nc3Error::verbose_nonfatal); - Nc3File dataFile(filename, Nc3File::ReadOnly); + + Nc3File dataFile(filename.c_str(), 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"))) + Nc3Var *lonVar, *latVar, *altVar, *timeVar; + std::cout << "before get_var(lon) " << std::endl; + if (!(lonVar = dataFile.get_var("longitude"))) return NC_ERR; - if (!(latVar = dataFile.get_var("lat"))) + std::cout << "before get_var(lat) " << std::endl; + if (!(latVar = dataFile.get_var("latitude"))) return NC_ERR; - if (!(altVar = dataFile.get_var("z"))) + std::cout << "before get_var(z) " << std::endl; + if (!(altVar = dataFile.get_var("altitude"))) return NC_ERR; + std::cout << "before get_var(time) " << std::endl; + if (!(timeVar= dataFile.get_var("time"))) + return NC_ERR; + std::cout << "after get_var(time) " << std::endl; // Get the lat/lon data from the file. if (!lonVar->get(longitude, NLON)) @@ -95,35 +103,36 @@ int NetCDF_XYZ::readNetCDF(const char* filename) { // 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"))) + if (!(uVar = dataFile.get_var("U"))) return NC_ERR; - if (!(vVar = dataFile.get_var("v"))) + if (!(vVar = dataFile.get_var("V"))) return NC_ERR; - if (!(wVar = dataFile.get_var("w"))) + if (!(wVar = dataFile.get_var("W"))) return NC_ERR; - if (!(dudxVar = dataFile.get_var("dudx"))) + if (!(dudxVar = dataFile.get_var("DUDX"))) return NC_ERR; - if (!(dvdxVar = dataFile.get_var("dvdx"))) + if (!(dvdxVar = dataFile.get_var("DVDX"))) return NC_ERR; - if (!(dwdxVar = dataFile.get_var("dwdx"))) + if (!(dwdxVar = dataFile.get_var("DWDX"))) return NC_ERR; - if (!(dudyVar = dataFile.get_var("dudy"))) + if (!(dudyVar = dataFile.get_var("DUDY"))) return NC_ERR; - if (!(dvdyVar = dataFile.get_var("dvdy"))) + if (!(dvdyVar = dataFile.get_var("DVDY"))) return NC_ERR; - if (!(dwdyVar = dataFile.get_var("dwdy"))) + if (!(dwdyVar = dataFile.get_var("DWDY"))) return NC_ERR; - if (!(dudzVar = dataFile.get_var("dudz"))) + if (!(dudzVar = dataFile.get_var("DUDZ"))) return NC_ERR; - if (!(dvdzVar = dataFile.get_var("dvdz"))) + if (!(dvdzVar = dataFile.get_var("DVDZ"))) return NC_ERR; - if (!(dwdzVar = dataFile.get_var("dwdz"))) + if (!(dwdzVar = dataFile.get_var("DWDZ"))) return NC_ERR; - if (!(trbVar = dataFile.get_var("trb"))) - return NC_ERR; - if (!(dpibdxVar = dataFile.get_var("dpibdx"))) + //std::cout << "before get_var(trb) " << std::endl; + if (!(trbVar = dataFile.get_var("TRB"))) + return NC_ERR; + if (!(dpibdxVar = dataFile.get_var("DPIBDX"))) return NC_ERR; - if (!(dpibdyVar = dataFile.get_var("dpibdy"))) + if (!(dpibdyVar = dataFile.get_var("DPIBDY"))) return NC_ERR; if (!uVar->set_cur(NREC, 0, 0, 0)) @@ -151,11 +160,11 @@ int NetCDF_XYZ::readNetCDF(const char* filename) { if (!dwdzVar->set_cur(NREC, 0, 0, 0)) return NC_ERR; if (!trbVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; + return NC_ERR; if (!dpibdxVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; + return NC_ERR; if (!dpibdyVar->set_cur(NREC, 0, 0, 0)) - return NC_ERR; + return NC_ERR; if (!uVar->get(u, 1, NALT, NLAT, NLON)) return NC_ERR; @@ -181,17 +190,31 @@ int NetCDF_XYZ::readNetCDF(const char* filename) { 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; - + 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 1; // return an error because it is missing the trb, dpidx, dpidy variables return 0; } - +int NetCDF_XYZ::getValue(const std::string& varName) +{ + int value_out; + std::cout << "Inside getValue" << std::endl; + if (varName == "time") { + value_out = obtime; + } else { + std::cout << "Requested Variable unknown. \n"; + std::cout << varName << "\n"; + return false; + } + return value_out; +} + double NetCDF_XYZ::getValue(const int &i,const int &j,const int &k, const std::string& varName) { //Returned values are all in SI units @@ -314,7 +337,7 @@ double NetCDF_XYZ::calc_C(const int &i,const int &j,const int &k) 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 diff --git a/src/NetCDF_XYZ.h b/src/NetCDF_XYZ.h index 2df6247..bad1bc7 100644 --- a/src/NetCDF_XYZ.h +++ b/src/NetCDF_XYZ.h @@ -24,8 +24,9 @@ class NetCDF_XYZ : public NetCDF ~NetCDF_XYZ(); - int readNetCDF(const char* filename); + int readNetCDF(std::string filename); double getValue(const int &i,const int &j,const int &k,const std::string& varName); + int getValue(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); @@ -39,6 +40,7 @@ class NetCDF_XYZ : public NetCDF float* longitude; float* latitude; float* altitude; + int obtime; float* u; float* v; float* w; diff --git a/src/Params.cc b/src/Params.cc index 6e6f1ff..7684557 100644 --- a/src/Params.cc +++ b/src/Params.cc @@ -741,15 +741,15 @@ tt->enum_def.fields[1].val = MODE_RTZ; tt->single_val.e = MODE_XYZ; tt++; - + // Parameter 'analysis_type' // ctype is '_analysis_type_t' memset(tt, 0, sizeof(TDRPtable)); tt->ptype = ENUM_TYPE; tt->param_name = tdrpStrDup("analysis_type"); - tt->descr = tdrpStrDup("Type of analysis to perform"); - tt->help = tdrpStrDup("TODO"); + tt->descr = tdrpStrDup("Analysis type"); + tt->help = tdrpStrDup("WIND for Samurai wind analysis, THERMO for Thermo analysis from samurai output file, WIND_THERMO for Wind and then Thermo analysis (in-memory)"); tt->val_offset = (char *) &analysis_type - &_start_; tt->enum_def.name = tdrpStrDup("analysis_type_t"); tt->enum_def.nfields = 3; @@ -763,7 +763,6 @@ tt->enum_def.fields[2].val = WIND_THERMO; tt->single_val.e = WIND; tt++; - // Parameter 'projection' // ctype is '_projection_t' @@ -2111,6 +2110,66 @@ tt->single_val.f = 2.5; tt++; + // Parameter 'thermo_A_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("thermo_A_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->val_offset = (char *) &thermo_A_error - &_start_; + tt->single_val.f = 0.01; + tt++; + + // Parameter 'thermo_B_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("thermo_B_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->val_offset = (char *) &thermo_B_error - &_start_; + tt->single_val.f = 0.01; + tt++; + + // Parameter 'thermo_C_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("thermo_C_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->val_offset = (char *) &thermo_C_error - &_start_; + tt->single_val.f = 0.01; + tt++; + + // Parameter 'thermo_D_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("thermo_D_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->val_offset = (char *) &thermo_D_error - &_start_; + tt->single_val.f = 0.01; + tt++; + + // Parameter 'thermo_E_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("thermo_E_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->val_offset = (char *) &thermo_E_error - &_start_; + tt->single_val.f = 0.01; + tt++; + // Parameter 'qscat_rhov_error' // ctype is 'float' @@ -2591,6 +2650,18 @@ tt->single_val.b = pFALSE; tt++; + // Parameter 'wind_file' + // ctype is 'char*' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = STRING_TYPE; + tt->param_name = tdrpStrDup("wind_file"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup("Need for running Samurai in THERMO mode"); + tt->val_offset = (char *) &wind_file - &_start_; + tt->single_val.s = tdrpStrDup(""); + tt++; + // Parameter 'Comment 14' memset(tt, 0, sizeof(TDRPtable)); @@ -2820,6 +2891,66 @@ tt->array_vals[1].f = 1; tt++; + // Parameter 'bg_pip_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("bg_pip_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->array_offset = (char *) &_bg_pip_error - &_start_; + tt->array_n_offset = (char *) &bg_pip_error_n - &_start_; + tt->is_array = TRUE; + tt->array_len_fixed = FALSE; + tt->array_elem_size = sizeof(float); + tt->array_n = 2; + tt->array_vals = (tdrpVal_t *) + tdrpMalloc(tt->array_n * sizeof(tdrpVal_t)); + tt->array_vals[0].f = 10; + tt->array_vals[1].f = 10; + tt++; + + // Parameter 'bg_thetarhop_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("bg_thetarhop_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->array_offset = (char *) &_bg_thetarhop_error - &_start_; + tt->array_n_offset = (char *) &bg_thetarhop_error_n - &_start_; + tt->is_array = TRUE; + tt->array_len_fixed = FALSE; + tt->array_elem_size = sizeof(float); + tt->array_n = 2; + tt->array_vals = (tdrpVal_t *) + tdrpMalloc(tt->array_n * sizeof(tdrpVal_t)); + tt->array_vals[0].f = 10; + tt->array_vals[1].f = 10; + tt++; + + // Parameter 'bg_ftheta_error' + // ctype is 'float' + + memset(tt, 0, sizeof(TDRPtable)); + tt->ptype = FLOAT_TYPE; + tt->param_name = tdrpStrDup("bg_ftheta_error"); + tt->descr = tdrpStrDup(""); + tt->help = tdrpStrDup(""); + tt->array_offset = (char *) &_bg_ftheta_error - &_start_; + tt->array_n_offset = (char *) &bg_ftheta_error_n - &_start_; + tt->is_array = TRUE; + tt->array_len_fixed = FALSE; + tt->array_elem_size = sizeof(float); + tt->array_n = 2; + tt->array_vals = (tdrpVal_t *) + tdrpMalloc(tt->array_n * sizeof(tdrpVal_t)); + tt->array_vals[0].f = 10; + tt->array_vals[1].f = 10; + tt++; + // Parameter 'i_filter_length' // ctype is 'float' @@ -2999,7 +3130,7 @@ tt->array_vals[0].f = -1; tt->array_vals[1].f = -1; tt++; - + // Parameter 'Comment 15' memset(tt, 0, sizeof(TDRPtable)); @@ -3008,208 +3139,223 @@ tt->comment_hdr = tdrpStrDup("VARIABLES NEEDED BY THERMO"); tt->comment_text = tdrpStrDup(""); tt++; - + // Parameter 'i_pip_bcL' - // ctype is 'int' + // ctype is 'char*' memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("i_pip_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &i_pip_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'i_pip_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("i_pip_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &i_pip_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'i_thetarhop_bcL' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("i_thetarhop_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &i_thetarhop_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'i_thetarhop_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("i_thetarhop_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &i_thetarhop_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'i_ftheta_bcL' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("i_ftheta_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &i_ftheta_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'i_ftheta_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("i_ftheta_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &i_ftheta_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'j_pip_bcL' - // ctype is 'int' + // ctype is 'char*' memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("j_pip_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &j_pip_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'j_pip_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("j_pip_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &j_pip_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'j_thetarhop_bcL' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("j_thetarhop_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &j_thetarhop_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'j_thetarhop_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("j_thetarhop_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &j_thetarhop_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'j_ftheta_bcL' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("j_ftheta_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &j_ftheta_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'j_ftheta_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("j_ftheta_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &j_ftheta_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'k_pip_bcL' - // ctype is 'int' + // ctype is 'char*' memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("k_pip_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &k_pip_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'k_pip_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("k_pip_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &k_pip_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'k_thetarhop_bcL' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("k_thetarhop_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &k_thetarhop_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'k_thetarhop_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("k_thetarhop_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &k_thetarhop_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'k_ftheta_bcL' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("k_ftheta_bcL"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &k_ftheta_bcL - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // Parameter 'k_ftheta_bcR' - // ctype is 'int' + // ctype is 'char*' + memset(tt, 0, sizeof(TDRPtable)); - tt->ptype = INT_TYPE; + tt->ptype = STRING_TYPE; tt->param_name = tdrpStrDup("k_ftheta_bcR"); tt->descr = tdrpStrDup(""); tt->help = tdrpStrDup(""); tt->val_offset = (char *) &k_ftheta_bcR - &_start_; - tt->single_val.i = -999; + tt->single_val.s = tdrpStrDup("R0"); tt++; - + // trailing entry has param_name set to NULL tt->param_name = NULL; diff --git a/src/Params.hh b/src/Params.hh index 14a725f..7405611 100644 --- a/src/Params.hh +++ b/src/Params.hh @@ -622,6 +622,16 @@ public: float qscat_rhou_error; + float thermo_A_error; + + float thermo_B_error; + + float thermo_C_error; + + float thermo_D_error; + + float thermo_E_error; + float qscat_rhov_error; float ascat_rhou_error; @@ -696,6 +706,8 @@ public: tdrp_bool_t use_fractl_errors; + char* wind_file; + float *_mc_weight; int mc_weight_n; @@ -729,6 +741,15 @@ public: float *_bg_tempk_error; int bg_tempk_error_n; + float *_bg_pip_error; + int bg_pip_error_n; + + float *_bg_thetarhop_error; + int bg_thetarhop_error_n; + + float *_bg_ftheta_error; + int bg_ftheta_error_n; + float *_i_filter_length; int i_filter_length_n; @@ -756,41 +777,41 @@ public: float *_k_max_wavenumber; int k_max_wavenumber_n; - int i_pip_bcL; + char* i_pip_bcL; + + char* i_pip_bcR; - int i_pip_bcR; + char* i_thetarhop_bcL; - int i_thetarhop_bcL; + char* i_thetarhop_bcR; - int i_thetarhop_bcR; + char* i_ftheta_bcL; - int i_ftheta_bcL; + char* i_ftheta_bcR; - int i_ftheta_bcR; - - int j_pip_bcL; + char* j_pip_bcL; - int j_pip_bcR; + char* j_pip_bcR; - int j_thetarhop_bcL; + char* j_thetarhop_bcL; - int j_thetarhop_bcR; + char* j_thetarhop_bcR; - int j_ftheta_bcL; + char* j_ftheta_bcL; - int j_ftheta_bcR; + char* j_ftheta_bcR; - int k_pip_bcL; + char* k_pip_bcL; - int k_pip_bcR; + char* k_pip_bcR; - int k_thetarhop_bcL; + char* k_thetarhop_bcL; - int k_thetarhop_bcR; + char* k_thetarhop_bcR; - int k_ftheta_bcL; + char* k_ftheta_bcL; - int k_ftheta_bcR; + char* k_ftheta_bcR; char _end_; // end of data region // needed for zeroing out data @@ -799,7 +820,7 @@ private: void _init(); - mutable TDRPtable _table[210]; + mutable TDRPtable _table[220]; const char *_className; diff --git a/src/VarDriver.cpp b/src/VarDriver.cpp index 1b311f6..cddbd02 100644 --- a/src/VarDriver.cpp +++ b/src/VarDriver.cpp @@ -1671,6 +1671,11 @@ bool VarDriver::parseXMLconfig(const XMLNode& config) configKeys.insert("sfmr_windspeed_error"); configKeys.insert("qscat_rhou_error"); configKeys.insert("qscat_rhov_error"); + configKeys.insert("thermo_A_error"); + configKeys.insert("thermo_B_error"); + configKeys.insert("thermo_C_error"); + configKeys.insert("thermo_D_error"); + configKeys.insert("thermo_E_error"); configKeys.insert("ascat_rhou_error"); configKeys.insert("ascat_rhov_error"); configKeys.insert("amv_rhou_error"); @@ -1688,6 +1693,9 @@ bool VarDriver::parseXMLconfig(const XMLNode& config) configKeys.insert("bg_qv_error"); configKeys.insert("bg_rhoa_error"); configKeys.insert("bg_qr_error"); + configKeys.insert("bg_pip_error"); + configKeys.insert("bg_thetarhop_error"); + configKeys.insert("bg_ftheta_error"); if (fixedGrid) { configKeys.insert("ref_time"); diff --git a/src/VarDriver.h b/src/VarDriver.h index f31529d..6cc6f71 100644 --- a/src/VarDriver.h +++ b/src/VarDriver.h @@ -84,6 +84,11 @@ class VarDriver //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); + enum analysModes { + WIND, + THERMO, + WIND_THERMO + }; protected: bool fixedGrid; // Indicates if the grid dims come from the config file or run call @@ -134,6 +139,7 @@ class VarDriver hdob, ict }; + Projection projection; bool read_met_obs_file(int suffix, std::string &filename, std::vector* metObVector); diff --git a/src/VarDriver3D.cpp b/src/VarDriver3D.cpp index b5f442f..138588a 100644 --- a/src/VarDriver3D.cpp +++ b/src/VarDriver3D.cpp @@ -50,7 +50,7 @@ VarDriver3D::~VarDriver3D() bool VarDriver3D::initialize() { // Run a 3D vortex background field - cout << "Initializing SAMURAI 3D" << endl; + std::cout << "Initializing SAMURAI 3D" << std::endl; return validateDriver(); } @@ -59,7 +59,7 @@ bool VarDriver3D::initialize() bool VarDriver3D::initialize(const XMLNode& configuration) { // Run a 3D vortex background field - cout << "Initializing SAMURAI 3D" << endl; + std::cout << "Initializing SAMURAI 3D" << std::endl; // Parse the XML configuration file if (!parseXMLconfig(configuration)) return false; @@ -72,7 +72,7 @@ bool VarDriver3D::initialize(const XMLNode& configuration) bool VarDriver3D::initialize(const samurai_config &configSam) { // Run a 3D vortex background field - cout << "Initializing SAMURAI 3D" << endl; + std::cout << "Initializing SAMURAI 3D" << std::endl; // Parse the Samurai configuration structure passed from COAMPS if (! parseSamuraiConfig(configSam)) return false; @@ -98,30 +98,25 @@ bool VarDriver3D::validateDriver() } else if (configHash["mode"] == "RTZ") { runMode = RTZ; } else { - cout << "Unrecognized run mode " << configHash["mode"] << ", Aborting...\n"; + std::cout << "Unrecognized run mode " << configHash["mode"] << ", Aborting...\n"; return false; } // Print the analysis_type from the TDRP config file - std::cout << "Analysis type: " << configHash["analysis_type"] << std::endl; - if (configHash["analysis_type"] == "WIND") - { - numVars = 7; // Number of variables on which to perform the anslysis - obMetaSize = 7; // Size of the observation Meta data - } - if (configHash["analysis_type"] == "THERMO") - { - numVars = 3; // Number of variables on which to perform the anslysis - obMetaSize = 7; // Size of the observation Meta data + //std::string analysis_type = configHash["analysis_type"]; + + if(configHash["analysis_type"] == "WIND") { + numVars = 7; // Number of variables on which to perform the anslysis + obMetaSize = 7; // Size of the observation Meta data + analysisMode = VarDriver::WIND; + } else if (configHash["analysis_type"] == "THERMO") { + numVars = 3; // Number of variables on which to perform the anslysis + obMetaSize = 7; // Size of the observation Meta data + analysisMode = VarDriver::THERMO; + } else { + std::cout << "Currently unsupported Analysis type: " << configHash["analysis_type"] << ", Aborting..." << std::endl; + return false; } - if(configHash["analysis_type"] != "WIND") - { - std::cout << "i_pip_bcL = " << configHash["i_pip_bcL"] << std::endl; - std::cout << "j_pip_bcR = " << configHash["j_pip_bcR"] << std::endl; - std::cout << "k_ftheta_bcL = " << configHash["k_ftheta_bcL"] << std::endl; - std::cout << "Currently unsupported Analysis type: " << configHash["analysis_type"] << ", Aborting..." << std::endl; - return false; - } bool fractlBkgd = ( configHash["bkgd_obs_interpolation"] == "fractl" ); bool loadBG = ( configHash["load_background"] == "true" ); @@ -149,7 +144,7 @@ bool VarDriver3D::validateDriver() // that's not fully implemented in some compilers yet, so for the moment we're just assuming it's a directory dataPath = configHash["data_directory"]; if (!DirectoryExists(dataPath)) { - std::cout << "Can't find data directory: " << configHash["data_directory"] << endl; + std::cout << "Can't find data directory: " << configHash["data_directory"] << std::endl; return false; } @@ -157,7 +152,7 @@ bool VarDriver3D::validateDriver() std::string outputPath(configHash["output_directory"]); if (!DirectoryExists(outputPath)) { - std::cout << "Can't find output directory: " << configHash["output_directory"] << endl; + std::cout << "Can't find output directory: " << configHash["output_directory"] << std::endl; return false; } @@ -209,17 +204,20 @@ bool VarDriver3D::gridDependentInit() std::string refSounding = configHash["ref_state"]; refstate = new ReferenceState(refSounding); - // cout << "Reference profile: Z\t\tQv\tRhoa\tRho\tH\tTemp\tPressure\n"; + + + + // std::cout << "Reference profile: Z\t\tQv\tRhoa\tRho\tH\tTemp\tPressure\n"; for (real k = kmin; k < kmax + kincr; k += kincr) { - // cout << " " << k << "\t"; + // std::cout << " " << k << "\t"; for (int i = 0; i < 6; i++) { real var = refstate->getReferenceVariable(i, k * 1000); if (i == 0) var = refstate->bhypInvTransform(var); - // cout << setw(9) << setprecision(4) << var << "\t"; + // std::cout << setw(9) << setprecision(4) << var << "\t"; } - // cout << "\n"; + // std::cout << "\n"; } - cout << setprecision(9); + std::cout << setprecision(9); // Set the maximum number of iterations to the multipass reduction factor // Multiple outer loops will reduce the cutoff wavelengths and background error variance @@ -235,13 +233,18 @@ bool VarDriver3D::gridDependentInit() if ( fixedGrid) { readFrameCenters(); if ( ! findReferenceCenter() ) { - cout << "Error finding reference time, please check date and time in XML file\n"; + std::cout << "Error finding reference time, please check date and time in XML file\n"; return false; } } - // These are used to process the obs (bkg and met) + if (analysisMode == VarDriver::THERMO) { + std::string windFile = configHash["wind_file"]; + this->loadObservations(windFile,idim,jdim,kdim); + } + // These are used to process the obs (bkg and met) + uStateSize = 8 * (idim + 1) * (jdim + 1) * (kdim + 1) * (numVars); // bgU size bStateSize = (idim + 2) * (jdim + 2) * (kdim + 2) * numVars; // 2 mish points between nodes @@ -268,6 +271,7 @@ bool VarDriver3D::gridDependentInit() // Optionally load a set of background estimates and interpolate to the Gaussian mish +if (analysisMode == WIND) { int numbgObs = 0; if(bkgdAdapter != NULL) { START_TIMER(timeb); @@ -275,7 +279,7 @@ bool VarDriver3D::gridDependentInit() PRINT_TIMER("loadBackgroundObs", timeb); if (numbgObs < 0) { - cout << "Error loading background Obs\n"; + std::cout << "Error loading background Obs\n"; return false; } //NCAR - cleanup from dangling bkgAdapter pointer, as it's not needed elsewhere: @@ -288,23 +292,55 @@ bool VarDriver3D::gridDependentInit() std::string adjustBG = configHash["adjust_background"]; if ((adjustBG == "true") and numbgObs) { - if ( ! adjustBackground()) { - cout << "Error adjusting background\n"; + if ( ! adjustBackground(analysisMode)) { + std::cout << "Error adjusting background\n"; return false; } } +} - START_TIMER(timem); - if ( ! loadMetObs() ) - return false; - PRINT_TIMER("loadMetObs", timem); + if (analysisMode == WIND) { + START_TIMER(timem); + if ( ! loadMetObs() ) return false; + PRINT_TIMER("loadMetObs", timem); + } else if (analysisMode == VarDriver::THERMO) { + // Load the observations into a vector + obs = new real[obVector.size()*(obMetaSize+numVars*numDerivatives)]; + for (int m=0; m < obVector.size(); m++) { + int n = m*(obMetaSize+numVars*numDerivatives); + Observation ob = obVector.at(m); + obs[n] = ob.getOb(); + real invError = ob.getInverseError(); + if (!invError) { + cout << "Undefined instrument error specification for " << ob.getType() << " instrument type!\n"; + return false; + } + obs[n+1] = invError; + if (runMode == XYZ) { + obs[n+2] = ob.getCartesianX(); + obs[n+3] = ob.getCartesianY(); + } else if (runMode == RTZ) { + obs[n+2] = ob.getRadius(); + obs[n+3] = ob.getTheta(); + } + obs[n+4] = ob.getAltitude(); + obs[n+5] = ob.getType(); + obs[n+6] = ob.getTime(); + for (unsigned int var = 0; var < numVars; var++) { + for (unsigned int d = 0; d < numDerivatives; ++d) { + int wgt_index = n + obMetaSize + numVars*d + var; + obs[wgt_index] = ob.getWeight(var, d); + + } + } + } + } START_TIMER(timec); initObCost3D(); PRINT_TIMER("initObCost3D", timec); PRINT_TIMER("gridDependentInit", timei); - return true; } @@ -323,7 +359,7 @@ bool VarDriver3D::findReferenceCenter() // Note - we can't insert, since it already exists.. so we're doing the awkward 'GetMap', then assigning // a new value directly. This can definitely be cleaner, but let's get it correct first, clean later. configHash.GetMap()->find("ref_time")->second = std::to_string(Date(frametime)); - cout << "Found matching reference time " << PrintTime(reftime) << " at " << frameVector[fi].getLat() << ", " << frameVector[fi].getLon() << "\n"; + std::cout << "Found matching reference time " << PrintTime(reftime) << " at " << frameVector[fi].getLat() << ", " << frameVector[fi].getLon() << "\n"; return true; } } @@ -409,7 +445,7 @@ bool VarDriver3D::run() } else { configHash.update("save_mish", "false"); } - cout << "Outer Loop Iteration: " << iter << endl; + std::cout << "Outer Loop Iteration: " << iter << endl; START_TIMER(timei); obCost3D->initState(iter); PRINT_TIMER("Cost3D Init", timei); @@ -451,6 +487,7 @@ bool VarDriver3D::preProcessMetObs() { GPTLstart("VarDriver3D::preprocessMetObs"); + int pSize=0; vector rhoP; // Convert the bg dBZ back to Z for further processing with real radar data @@ -503,7 +540,7 @@ bool VarDriver3D::preProcessMetObs() iter++; } zeroClevel = height; - cout << "Found zero C level at " << zeroClevel << " based on reference sounding" << endl; + std::cout << "Found zero C level at " << zeroClevel << " based on reference sounding" << std::endl; // Load Met. Observations auto filenames = FileList(dataPath); @@ -511,7 +548,7 @@ bool VarDriver3D::preProcessMetObs() int processedFiles = 0; int attemptedFiles = 0; std::vector* metData = new std::vector; - cout << "Found " << filenames.size() << " data files to read..." << endl; + std::cout << "Found " << filenames.size() << " data files to read..." << std::endl; int totalFiles = filenames.size(); for (std::size_t i = 0; i < filenames.size(); ++i) { @@ -533,12 +570,13 @@ bool VarDriver3D::preProcessMetObs() suffix = "cfrad"; } - cout << "Processing " << file << " of type " << suffix << endl; + std::cout << "Processing " << file << " of type " << suffix << std::endl; attemptedFiles++; // Read different types of files std::string fullpath = dataPath + "/" + file; if (! read_met_obs_file(dataSuffix[suffix], fullpath, metData)) continue; + // std::cout << "Number of potential observations: " << metData->size() << std::endl; processedFiles++; @@ -560,7 +598,7 @@ bool VarDriver3D::preProcessMetObs() datetime endTime_ob = frameVector.back().getTime(); auto endTime = Date(endTime_ob); int prevobs = obVector.size(); - std::cout << "i = " << i << endl; + std::cout << "i = " << i << std::endl; for (std::size_t i = 0; i < metData->size(); ++i) { // Make sure the ob is within the time limits @@ -579,11 +617,11 @@ bool VarDriver3D::preProcessMetObs() // std::cout << "tcstart: " << PrintDate(startTime_ob) << ", tcend: " << PrintDate(endTime_ob) << ", obTime: " << PrintDate(obTime_ob) << std::endl; continue; } - // std::cout << "timeProblem = " << timeProblem << endl; + // std::cout << "timeProblem = " << timeProblem << std::endl; int fi = std::chrono::duration_cast(obTime_ob - startTime_ob).count(); - // std::cout << "fi = " << fi << endl; + // std::cout << "fi = " << fi << std::endl; // if ((fi < 0) or (fi > (int)frameVector.size())) { - // cout << "**Time problem with observation " << fi << ", " << startTime << ", " << obTime << endl; + // std::cout << "**Time problem with observation " << fi << ", " << startTime << ", " << obTime << std::endl; // timeProblem++; // continue; // } testing @@ -676,7 +714,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = rho*((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["dropsonde_rhou_error"])); obVector.push_back(varOb); @@ -825,7 +863,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = rho*((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["insitu_rhou_error"])); obVector.push_back(varOb); @@ -929,7 +967,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = ((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["qscat_rhou_error"])); obVector.push_back(varOb); @@ -962,7 +1000,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = ((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["ascat_rhou_error"])); obVector.push_back(varOb); @@ -995,7 +1033,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = ((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["amv_rhou_error"])); obVector.push_back(varOb); @@ -1321,7 +1359,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = rho*((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["mesonet_rhou_error"])); obVector.push_back(varOb); @@ -1398,7 +1436,7 @@ bool VarDriver3D::preProcessMetObs() } else if (runMode == RTZ) { rhou = rho*((u - Um)*obX + (v - Vm)*obY)/obRadius; } - //cout << "RhoU: " << rhou << endl; + //std::cout << "RhoU: " << rhou << std::endl; varOb.setOb(rhou); varOb.setError(std::stof(configHash["aeri_rhou_error"])); obVector.push_back(varOb); @@ -1826,20 +1864,22 @@ bool VarDriver3D::preProcessMetObs() // Show a summary of what got tossed out - cout << "Observation problem: " << obsProblem << ", Time problem: " << timeProblem + std::cout << "Observation problem: " << obsProblem << ", Time problem: " << timeProblem << ", Coordinate problem: " << coordProblem << ", Domain problem: " << domainProblem - << ", Radius problem: " << radiusProblem << endl; + << ", Radius problem: " << radiusProblem << std::endl; + std::cout << "Potential size: " << pSize << std::endl; std::cout << "obVector size: " << obVector.size() << std::endl; + pSize += (metData->size() - (obsProblem+timeProblem+coordProblem+domainProblem+radiusProblem)); int newobs = obVector.size() - prevobs; // if (metData->size() > 0) { if (newobs > 0) { - cout << "Processed " << newobs << " observations from " << metData->size() << " entries (" + std::cout << "Processed " << newobs << " observations from " << metData->size() << " entries (" << 100.0*(float)newobs/(6*(float)metData->size()) << "%) file: " << file << std::endl; } else { - cout << "No valid observations in file\n"; + std::cout << "No valid observations in file\n"; } - cout << obVector.size() << " total observations." << " ( " << attemptedFiles << " of " << totalFiles << " files processed ) " << endl; + std::cout << obVector.size() << " total observations." << " ( " << attemptedFiles << " of " << totalFiles << " files processed ) " << std::endl; } delete metData; @@ -1968,7 +2008,7 @@ bool VarDriver3D::preProcessMetObs() } } - cout << obVector.size() << " total observations including pseudo-obs for W and mass continuity" << endl; + std::cout << obVector.size() << " total observations including pseudo-obs for W and mass continuity" << std::endl; #if IO_WRITEOBS GPTLstart("VarDriver3D::preprocessMetObs->writeobs"); @@ -1989,7 +2029,7 @@ bool VarDriver3D::preProcessMetObs() *os++ = "Weight 4"; *os++ = "Weight 5"; *os++ = "Weight 6"; - obstream << endl; */ + obstream << std::endl; */ ostream_iterator od(obstream, "\t "); ostream_iterator oi(obstream, "\t "); @@ -1998,7 +2038,7 @@ bool VarDriver3D::preProcessMetObs() *od++ = ob.getOb(); real invError = ob.getInverseError(); if (!invError) { - cout << "Undefined instrument error specification for " << ob.getType() << "instrument type!\n"; + std::cout << "Undefined instrument error specification for " << ob.getType() << "instrument type!\n"; return false; } *od++ = invError; @@ -2017,7 +2057,7 @@ bool VarDriver3D::preProcessMetObs() *od++ = ob.getWeight(var, d); } } - obstream << endl; + obstream << std::endl; } GPTLstop("VarDriver3D::preprocessMetObs->writeobs"); #endif @@ -2031,7 +2071,7 @@ bool VarDriver3D::preProcessMetObs() obs[n] = ob.getOb(); real invError = ob.getInverseError(); if (!invError) { - cout << "Undefined instrument error specification for " << ob.getType() << "instrument type!\n"; + std::cout << "Undefined instrument error specification for " << ob.getType() << "instrument type!\n"; return false; } obs[n+1] = invError; @@ -2055,16 +2095,211 @@ bool VarDriver3D::preProcessMetObs() // All done preprocessing if (!processedFiles) { - cout << "No files processed, nothing to do :(" << endl; + std::cout << "No files processed, nothing to do :(" << std::endl; // return 0; } else { - cout << "Finished preprocessing " << processedFiles << " files." << endl; + std::cout << "Finished preprocessing " << processedFiles << " files." << std::endl; } GPTLstop("VarDriver3D::preprocessMetObs"); return true; } + +//bool VarDriverThermo::loadObservations(std::string filename,const int &metFile_idim,const int &metFile_jdim,const int &metFile_kdim , QList* obVector) +bool VarDriver3D::loadObservations(std::string filename,const int &metFile_idim,const int &metFile_jdim,const int &metFile_kdim) +{ + + std::cout << "Reading in WIND output file: " << filename << std::endl; + if (runMode == XYZ) { + ncFile = new NetCDF_XYZ(metFile_idim,metFile_jdim,metFile_kdim); + } else if (runMode == RTZ) { + // ncFile = new NetCDF_RTZ(metFile_idim,metFile_jdim,metFile_kdim); + std::cout << "Runmode RTZ is not supported in VarDriver3D::loadObservations" << std::endl; + } + + std::cout << "Read in NetCDF File " << filename.c_str() << "..." << std::endl; + if (ncFile->readNetCDF(filename) != 0) { + std::cout << "Error reading NetCDF file\n"; + exit(1); + } + + std::cout << "Load Observations ... " << std::endl; + + int nalt = metFile_kdim; + int nx = metFile_idim; + int ny = metFile_jdim; + + //datetime obTime; + + // NCAR note: I'm still confused about times - sometimes they're reals, sometimes ints?? Treating as seconds here + // bgTime.setTime_t(time); + // bgTime.setTimeSpec(Qt::UTC); + //std::cout << "iTime: " << itime << std::endl; + int64_t obTime = ncFile->getValue("time"); + //obTime = std::chrono::time_point(std::chrono::duration(foo))); + //obTime = std::chrono::time_point(std::chrono::duration(foo)); + + + std::cout << "loadObservations: timeVar: " << obTime << std::endl; + int fi = 1; + // std::cout << "loadObservations: timePoint: " << obTime.c_str() << std::endl; +#if 0 + QString file,datestr,timestr; + file = metFile.section("/",-1); + datestr = file.left(8); + QDate date = QDate::fromString(datestr, "yyyyMMdd"); + timestr = file.section("_",-1).section(".",0,0); + QTime time; + if (timestr.size()==2) { + time = QTime::fromString(timestr, "HH"); + } else if (timestr.size()==4) { + time = QTime::fromString(timestr,"HHmm"); + } else { + std::cout << "Implement reading routine for filenames which don't look like yyyymmdd_hh.nc \n"; + exit(1); + } + + QDateTime obTime; + obTime = QDateTime(date, time, Qt::UTC); + QString obstring = obTime.toString(Qt::ISODate); + QDateTime startTime = frameVector.front().getTime(); + QDateTime endTime = frameVector.back().getTime(); + QString tcstart = startTime.toString(Qt::ISODate); + QString tcend = endTime.toString(Qt::ISODate); + int fi = startTime.secsTo(obTime); +#endif + if ((fi < 0) or (fi > (int)frameVector.size())) { + std::cout << "Time problem with observation " << "size(frameVector): " << (int)frameVector.size() << " " << fi << std::endl; + exit(1); + } + + + //if (configHash["mode"] == "RTZ") { + //} else if (configHash["mode"] == "XYZ") { + if (runMode == RTZ) { + std::cout << "RTZ run mode not implemented yet " << std::endl; + + } else if (runMode == XYZ) { + std::cout << "Entered XYZ run mode " << std::endl; + + int nlon = metFile_idim; + int nlat = metFile_jdim; + int nalt = metFile_kdim; + std::cout << "loadObservations: nlon*nlat*nalt " << nlon*nlat*nalt << std::endl; + for (int i = 0; i < nlon; ++i) { + for (int j = 0; j < nlat; ++j) { + for (int k = 0; k < nalt; ++k) { + + Observation varOb; + + double lon = ncFile->getValue(i,j,k,"lon"); + double lat = ncFile->getValue(i,j,k,"lat"); + double alt = ncFile->getValue(i,j,k,"z"); + double alt_km = alt/1000.0; + + //Checking if obs in domain is still missing! + // if ((r_km < imin) or (r_km > imax) or + // (lambda < jmin) or (lambda > jmax) or + // (alt_km < kmin) or (alt_km > kmax)) + // continue; + + // Geographic functions + GeographicLib::TransverseMercatorExact tm = GeographicLib::TransverseMercatorExact::UTM(); + double referenceLon = std::stof(configHash["ref_lon"]); + double tcX, tcY, cartX, cartY; + tm.Forward(referenceLon, frameVector[fi].getLat() , frameVector[fi].getLon() , tcX, tcY); + tm.Forward(referenceLon,lat,lon,cartX,cartY); + + real obX = (cartX - tcX)/1000.; + real obY = (cartY - tcY)/1000.; + varOb.setType(MetObs::model); + varOb.setCartesianX(obX); + varOb.setCartesianY(obY); + varOb.setAltitude(alt_km); + varOb.setTime(obTime); + + // Initialize the weights + for (unsigned int var = 0; var < numVars; ++var) { + for (unsigned int d = 0; d < numDerivatives; ++d) { + varOb.setWeight(0.0, var, d); + } + } + + double a,b,c,d,e; + + a = ncFile->calc_A(i,j,k); + b = ncFile->calc_B(i,j,k); + c = ncFile->calc_C(i,j,k); + d = ncFile->calc_D(i,j,k); + e = ncFile->calc_E(i,j,k); + + + double thetarhobar = ncFile->getValue(i,j,k,"trb"); + + double u = ncFile->getValue(i,j,k,"u"); + double v = ncFile->getValue(i,j,k,"v"); + double w = ncFile->getValue(i,j,k,"w"); + double wspd = u*u+v*v; + + + if (a==-999 or b==-999 or c==-999 or d==-999 or e==-999 or thetarhobar==-999){ + std::cout << "loadObservations: a,b,c,d,e " << a << " " << b << " " << c << " " << " " << d << " " << e << std::endl; + continue; + } + + if (a==-999 or b==-999 or c==-999 or d==-999 or e==-999 or thetarhobar==-999){std::cout << "Skip this ... \n";} + + float c_p = 1005.7; + float g = 9.81; + + varOb.setOb(a*1E3*1E3); + varOb.setWeight(-0.001*1E3,0,1); + varOb.setError(std::stof(configHash["thermo_A_error"])); + obVector.push_back(varOb); + varOb.setWeight(0,0,1); + + varOb.setOb(b*1E3*1E3); + varOb.setWeight(-0.001*1E3,0,2); + varOb.setError(std::stof(configHash["thermo_B_error"])); + obVector.push_back(varOb); + varOb.setWeight(0,0,2); + + varOb.setOb(c*1E3*1E3); + varOb.setWeight(-0.001*1E3,0,3); + varOb.setWeight(g*1E3*1E3/(c_p*thetarhobar*thetarhobar),1,0); + varOb.setError(std::stof(configHash["thermo_C_error"])); + obVector.push_back(varOb); + varOb.setWeight(0,0,3); + varOb.setWeight(0,1,0); + + varOb.setOb(d); + varOb.setWeight(1,1,1); + varOb.setError(std::stof(configHash["thermo_D_error"])); + obVector.push_back(varOb); + varOb.setWeight(0,1,1); + + varOb.setOb(e); + varOb.setWeight(1,1,2); + varOb.setError(std::stof(configHash["thermo_E_error"])); + obVector.push_back(varOb); + varOb.setWeight(0,1,2); + + } + } + } + + } else { + std::cout << "Unrecognized run mode " << configHash["mode"] << ", Aborting...\n"; + return false; + } + std::cout << "loadObservations: obVector.size(): " << obVector.size() << std::endl; + + return true; + +} + + /* Load the meteorological observations from a file into a vector */ bool VarDriver3D::loadPreProcessMetObs() @@ -2075,7 +2310,7 @@ bool VarDriver3D::loadPreProcessMetObs() real iPos, jPos, kPos, ob, error; int type; int64_t time; - cout << "Loading preprocessed observations from samurai_Observations.in" << endl; + std::cout << "Loading preprocessed observations from samurai_Observations.in" << std::endl; // Open and read the file std::string obFilename = dataPath + "/samurai_Observations.in"; @@ -2117,7 +2352,7 @@ bool VarDriver3D::loadPreProcessMetObs() obs[n] = ob.getOb(); real invError = ob.getInverseError(); if (!invError) { - cout << "Undefined instrument error specification for " << ob.getType() << "instrument type!\n"; + std::cout << "Undefined instrument error specification for " << ob.getType() << "instrument type!\n"; return false; } obs[n+1] = invError; @@ -2187,7 +2422,7 @@ int VarDriver3D::loadBackgroundObs() return bgIn.size() * obMetaSize / 11; } -bool VarDriver3D::adjustBackground() +bool VarDriver3D::adjustBackground(int analysisMode) { // Set the minimum filter length to the background resolution, not the analysis resolution // to avoid artifacts when running interpolating to small mesoscale grids @@ -2301,30 +2536,42 @@ bool VarDriver3D::adjustBackground() } // Store and set the background errors - std::string bgError[7]; - bgError[0] = configHash["bg_rhou_error"]; - bgError[1] = configHash["bg_rhov_error"]; - bgError[2] = configHash["bg_rhow_error"]; - bgError[3] = configHash["bg_tempk_error"]; - bgError[4] = configHash["bg_qv_error"]; - bgError[5] = configHash["bg_rhoa_error"]; - bgError[6] = configHash["bg_qr_error"]; + std::string bgError[numVars]; + if(analysisMode == WIND) { + bgError[0] = configHash["bg_rhou_error"]; + bgError[1] = configHash["bg_rhov_error"]; + bgError[2] = configHash["bg_rhow_error"]; + bgError[3] = configHash["bg_tempk_error"]; + bgError[4] = configHash["bg_qv_error"]; + bgError[5] = configHash["bg_rhoa_error"]; + bgError[6] = configHash["bg_qr_error"]; + } else if (analysisMode == THERMO) { + bgError[0] = configHash["bg_pip_error"]; + bgError[1] = configHash["bg_thetarhop_error"]; + bgError[2] = configHash["bg_ftheta_error"]; + } std::string bg_interpolation_error = "1.0"; if (configHash.exists("bg_interpolation_error")) { bg_interpolation_error = configHash["bg_interpolation_error"]; - cout << "Setting background interpolation error to " << bg_interpolation_error << "\n"; + std::cout << "Setting background interpolation error to " << bg_interpolation_error << "\n"; } else { - cout << "Using default background interpolation error of 1.0\n"; + std::cout << "Using default background interpolation error of 1.0\n"; + } + if(analysisMode == WIND) { + configHash.update("bg_rhou_error", bg_interpolation_error); + configHash.update("bg_rhov_error", bg_interpolation_error); + configHash.update("bg_rhow_error", bg_interpolation_error); + configHash.update("bg_tempk_error", bg_interpolation_error); + configHash.update("bg_qv_error", bg_interpolation_error); + configHash.update("bg_rhoa_error", bg_interpolation_error); + configHash.update("bg_qr_error", bg_interpolation_error); + } else if (analysisMode == THERMO) { + configHash.update("bg_pip_error", bg_interpolation_error); + configHash.update("bg_thetarhop_error", bg_interpolation_error); + configHash.update("bg_ftheta_error", bg_interpolation_error); } - configHash.update("bg_rhou_error", bg_interpolation_error); - configHash.update("bg_rhov_error", bg_interpolation_error); - configHash.update("bg_rhow_error", bg_interpolation_error); - configHash.update("bg_tempk_error", bg_interpolation_error); - configHash.update("bg_qv_error", bg_interpolation_error); - configHash.update("bg_rhoa_error", bg_interpolation_error); - configHash.update("bg_qr_error", bg_interpolation_error); - configHash.update("save_mish", "true"); + configHash.update("save_mish", "true"); // Adjust the background field to the spline mish @@ -2337,7 +2584,7 @@ bool VarDriver3D::adjustBackground() } else if (runMode == RTZ) { bgCost3D = new CostFunctionRTZ(projection, numbgObs, bStateSize); } - bgCost3D->initialize(&configHash, bgU, bgObs, refstate, numVars, numDerivatives, obMetaSize); + bgCost3D->initialize(&configHash, bgU, bgObs, refstate, numVars, numDerivatives, obMetaSize, analysisMode); // Set the iteration to zero -- // this will prevent writing the background file until after the adjustment @@ -2356,14 +2603,20 @@ bool VarDriver3D::adjustBackground() // Reset the background errors - configHash.update("bg_rhou_error", bgError[0]); - configHash.update("bg_rhov_error", bgError[1]); - configHash.update("bg_rhow_error", bgError[2]); - configHash.update("bg_tempk_error", bgError[3]); - configHash.update("bg_qv_error", bgError[4]); - configHash.update("bg_rhoa_error", bgError[5]); - configHash.update("bg_qr_error", bgError[6]); - configHash.update("save_mish", "false"); + if(analysisMode == WIND) { + configHash.update("bg_rhou_error", bgError[0]); + configHash.update("bg_rhov_error", bgError[1]); + configHash.update("bg_rhow_error", bgError[2]); + configHash.update("bg_tempk_error", bgError[3]); + configHash.update("bg_qv_error", bgError[4]); + configHash.update("bg_rhoa_error", bgError[5]); + configHash.update("bg_qr_error", bgError[6]); + configHash.update("save_mish", "false"); + } else if (analysisMode == THERMO) { + configHash.update("bg_pip_error", bgError[0]); + configHash.update("bg_thetarhop_error", bgError[1]); + configHash.update("bg_ftheta_error", bgError[2]); + } // Convert the dBZ back to Z for further processing @@ -2428,6 +2681,18 @@ void VarDriver3D::updateAnalysisParams(const int& iteration) val = configHash[key]; configHash.update("bg_qr_error", val); + key = "bg_pip_error_" + iter; + val = configHash[key]; + configHash.update("bg_pip_error", val); + + key = "bg_thetarhop_error_" + iter; + val = configHash[key]; + configHash.update("bg_thetarhop_error", val); + + key = "bg_ftheta_error_" + iter; + val = configHash[key]; + configHash.update("bg_ftheta_error", val); + key = "mc_weight_" + iter; val = configHash[key]; configHash.update("mc_weight", val); @@ -2556,7 +2821,7 @@ bool VarDriver3D::validateConfig() bool VarDriver3D::loadBackgroundCoeffs() { // Load a set of coefficients directly for the same grid - cout << "Loading previous coefficients from samurai_Coefficients.in" << endl; + std::cout << "Loading previous coefficients from samurai_Coefficients.in" << std::endl; std::ifstream bgFile(dataPath + "/samurai_Coefficients.in"); if (!bgFile.is_open()) @@ -2585,9 +2850,9 @@ bool VarDriver3D::loadBackgroundCoeffs() numbgCoeffs++; } } - cout << numbgCoeffs << " background coeffients loaded" << endl; + std::cout << numbgCoeffs << " background coeffients loaded" << std::endl; if (numbgCoeffs != bStateSize) { - cout << "Error loading background coefficients" << endl; + std::cout << "Error loading background coefficients" << std::endl; return false; } return true; @@ -2834,24 +3099,24 @@ bool VarDriver3D::validateGrid() // so less than 4 gridpoints will cause a memory fault if (idim < 4) { - cout << "i dimension is less than 4 gridpoints and recursive filter will fail. Aborting...\n"; + std::cout << "i dimension is less than 4 gridpoints and recursive filter will fail. Aborting...\n"; return false; } if (jdim < 4) { - cout << "j dimension is less than 4 gridpoints and recursive filter will fail. Aborting...\n"; + std::cout << "j dimension is less than 4 gridpoints and recursive filter will fail. Aborting...\n"; return false; } if (kdim < 4) { - cout << "k dimension is less than 4 gridpoints and recursive filter will fail. Aborting...\n"; + std::cout << "k dimension is less than 4 gridpoints and recursive filter will fail. Aborting...\n"; return false; } // Define the sizes of the arrays we are passing to the cost function - cout << "iMin\tiMax\tiIncr\tjMin\tjMax\tjIncr\tkMin\tkMax\tkIncr\n"; - cout << imin << "\t" << imax << "\t" << iincr << "\t"; - cout << jmin << "\t" << jmax << "\t" << jincr << "\t"; - cout << kmin << "\t" << kmax << "\t" << kincr << "\n\n"; + std::cout << "iMin\tiMax\tiIncr\tjMin\tjMax\tjIncr\tkMin\tkMax\tkIncr\n"; + std::cout << imin << "\t" << imax << "\t" << iincr << "\t"; + std::cout << jmin << "\t" << jmax << "\t" << jincr << "\t"; + std::cout << kmin << "\t" << kmax << "\t" << kincr << "\n\n"; return true; } @@ -2868,23 +3133,23 @@ bool VarDriver3D::loadMetObs() delete[] bgWeights; if (! success) { - cout << "Error pre-processing observations\n"; + std::cout << "Error pre-processing observations\n"; return false; } } else { if (!loadPreProcessMetObs()) { - cout << "Error loading observations\n"; + std::cout << "Error loading observations\n"; return false; } } if (obVector.size() == 0) { // No observations so quit - cout << "No observations loaded, unable to perform analysis.\n"; + std::cout << "No observations loaded, unable to perform analysis.\n"; return false; } else { - cout << "Number of New Observations: " << obVector.size() << endl; + std::cout << "Number of New Observations: " << obVector.size() << std::endl; } return true; } @@ -2892,7 +3157,7 @@ bool VarDriver3D::loadMetObs() bool VarDriver3D::initObCost3D() { if (runMode == XYZ) { - if (std::stof(configHash["output_pressure_increment"]) > 0) { + if (std::stof(configHash["output_pressure_increment"]) > 0) { obCost3D = new CostFunctionXYP(projection, obVector.size(), bStateSize); } else if (configHash["output_COAMPS"] == "true") { CostFunctionCOAMPS *cf = new CostFunctionCOAMPS(projection, obVector.size(), bStateSize); @@ -2905,7 +3170,7 @@ bool VarDriver3D::initObCost3D() obCost3D = new CostFunctionRTZ(projection, obVector.size(), bStateSize); } - obCost3D->initialize(&configHash, bgU, obs, refstate, numVars, numDerivatives, obMetaSize); + obCost3D->initialize(&configHash, bgU, obs, refstate, numVars, numDerivatives, obMetaSize, analysisMode); return true; } diff --git a/src/VarDriver3D.h b/src/VarDriver3D.h index 143d2a6..8252e7e 100644 --- a/src/VarDriver3D.h +++ b/src/VarDriver3D.h @@ -19,6 +19,8 @@ #include "CostFunctionCOAMPS.h" #include "MetObs.h" #include "FrameCenter.h" +#include "NetCDF.h" +#include "NetCDF_XYZ.h" #include "BkgdAdapter.h" #include "Xml.h" #include @@ -77,6 +79,7 @@ class VarDriver3D : public VarDriver private: + NetCDF* ncFile; typedef BSplineBase SplineBase; typedef BSpline SplineD; // This is also declared in ReferenceState.h @@ -100,6 +103,7 @@ class VarDriver3D : public VarDriver bool loadPreProcessMetObs(); bool loadBGfromFile(); bool loadBackgroundCoeffs(); + bool loadObservations(std::string filename, const int &metFile_idim, const int &metFile_jdim, const int &metFile_kdim); int loadBackgroundObs(const char *background_fname); int loadBackgroundObs(int nx, int ny, int nsigma, char *ctdg, int delta, int iter, // time elements @@ -113,7 +117,7 @@ class VarDriver3D : public VarDriver float *th1, float *p1); int loadBackgroundObs(); - bool adjustBackground(); + bool adjustBackground(int analysis_mode); // bool adjustBackground(const int& bStateSize); void updateAnalysisParams(const int& iteration); @@ -141,6 +145,7 @@ class VarDriver3D : public VarDriver int jdim; int kdim; int runMode; + int analysisMode; // These 2 control sizes of data structure. Pulling them here // since they were computed the same way in 2 different places. diff --git a/src/paramdef.samurai b/src/paramdef.samurai index 32f3929..40500be 100644 --- a/src/paramdef.samurai +++ b/src/paramdef.samurai @@ -111,7 +111,7 @@ typedef enum { paramdef enum analysis_type_t { p_default = WIND; p_descr = "Analysis type"; - p_help = "WIND for Samurai wind analysis, THERMO for Thermo analysis from samurai output file, WIND_THERMO for Wind and then Thermo analysis (in-memory)" + p_help = "WIND for Samurai wind analysis, THERMO for Thermo analysis from samurai output file, WIND_THERMO for Wind and then Thermo analysis (in-memory)"; } analysis_type; typedef enum { @@ -762,6 +762,36 @@ paramdef float { p_help = ""; } qscat_rhou_error; +paramdef float { + p_default = 0.01; + p_descr = ""; + p_help = ""; +} thermo_A_error; + +paramdef float { + p_default = 0.01; + p_descr = ""; + p_help = ""; +} thermo_B_error; + +paramdef float { + p_default = 0.01; + p_descr = ""; + p_help = ""; +} thermo_C_error; + +paramdef float { + p_default = 0.01; + p_descr = ""; + p_help = ""; +} thermo_D_error; + +paramdef float { + p_default = 0.01; + p_descr = ""; + p_help = ""; +} thermo_E_error; + paramdef float { p_default = 2.5; p_descr = ""; @@ -1000,12 +1030,19 @@ paramdef string { p_help = "Need bkgd_obs_interpolation set to INTERP_FRACTL to be used"; } fractl_nc_file; + paramdef boolean { p_default = false; p_desc = "Use Fractl generated errors instead of ..."; p_help = "Need bkgd_obs_interpolation set to 'fractl' to be used"; } use_fractl_errors; +paramdef string { + p_default = ""; + p_desc = "Output generated by Samurai in WIND mode"; + p_help = "Need for running Samurai in THERMO mode"; +} wind_file; + commentdef { p_header = "ITERATION DEPENDENT SECTION"; p_help = "All of these need as many entries as num_iterations"; @@ -1077,6 +1114,24 @@ paramdef float { p_help = ""; } bg_tempk_error[]; +paramdef float { + p_default = {10.0, 10.0}; + p_descr = ""; + p_help = ""; +} bg_pip_error[]; + +paramdef float { + p_default = {10.0, 10.0}; + p_descr = ""; + p_help = ""; +} bg_thetarhop_error[]; + +paramdef float { + p_default = {10.0, 10.0}; + p_descr = ""; + p_help = ""; +} bg_ftheta_error[]; + paramdef float { p_default = {4.0, 2.0}; p_descr = ""; @@ -1136,110 +1191,110 @@ commentdef { p_help = ""; } -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } i_pip_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } i_pip_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } i_thetarhop_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } i_thetarhop_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } i_ftheta_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } i_ftheta_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } j_pip_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } j_pip_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } j_thetarhop_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } j_thetarhop_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } j_ftheta_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } j_ftheta_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } k_pip_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } k_pip_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } k_thetarhop_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } k_thetarhop_bcR; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; } k_ftheta_bcL; -paramdef int { - p_default = -999; +paramdef string { + p_default = "R0"; p_descr = ""; p_help = ""; -} k_ftheta_bcR; \ No newline at end of file +} k_ftheta_bcR; diff --git a/src/samurai.h b/src/samurai.h index 472c59b..01fd5ce 100644 --- a/src/samurai.h +++ b/src/samurai.h @@ -59,6 +59,9 @@ struct samurai_config { float bg_qv_error; float bg_rhoa_error; float bg_qr_error; + float bg_pip_error; + float bg_thetarhop_error; + float bg_ftheta_error; float mc_weight; float neumann_u_weight; float neumann_v_weight; @@ -95,6 +98,11 @@ struct samurai_config { float sfmr_windspeed_error; float qscat_rhou_error; float qscat_rhov_error; + float thermo_A_error; + float thermo_B_error; + float thermo_C_error; + float thermo_D_error; + float thermo_E_error; float ascat_rhou_error; float ascat_rhov_error; float amv_rhou_error;