diff --git a/ncar_scripts/TDRP/beltrami.tdrp b/ncar_scripts/TDRP/beltrami.tdrp index e259e61..05b5b20 100644 --- a/ncar_scripts/TDRP/beltrami.tdrp +++ b/ncar_scripts/TDRP/beltrami.tdrp @@ -98,6 +98,21 @@ bkgd_obs_interpolation = INTERP_SPLINE; mode = MODE_XYZ; +///////////// mode //////////////////////////////////// +// Type of analysis. +// +// WIND for wind analysis +// THERMO for Thermo analysis using samurai output from a file +// WIND_THERMO for wind analysis and then Thermo analysis +// +// Type: enum +// Options: +// WIND +// THERMO +// WIND_THERMO + +analysis_type = WIND; + ///////////// projection ////////////////////////////// // // Projection. diff --git a/ncar_scripts/TDRP/beltramiTR.tdrp b/ncar_scripts/TDRP/beltramiTR.tdrp new file mode 100644 index 0000000..ac73bd8 --- /dev/null +++ b/ncar_scripts/TDRP/beltramiTR.tdrp @@ -0,0 +1,1237 @@ + +// Overwrite bgU with the content of this file. +// +// If set to a valid file, the bgU array will be overwritten. +// +// Type: string + +debug_bgu_overwrite = ""; + +///////////// debug_kd //////////////////////////////// +// +// Debug the KD tree. +// +// If set to a non-zero value val, dump val/step KD tree lookup. +// +// Type: int + +debug_kd = 0; + +///////////// debug_kd_step /////////////////////////// +// +// Step for decrementing debug_kd. +// +// debug_kd is decremented by this value. +// +// Type: int + +debug_kd_step = 0; + +//====================================================================== +// +// BACKGROUND SECTION. +// +//====================================================================== + +///////////// load_background ///////////////////////// +// +// Tell Samurai to load background observations. +// +// TODO. +// +// Type: boolean + +load_background = FALSE; + +///////////// load_bg_coefficients //////////////////// +// +// Tell Samurai to load background coefficients. +// +// TODO. +// +// Type: boolean + +load_bg_coefficients = FALSE; + +///////////// adjust_background /////////////////////// +// +// Tell Samurai to adjust the background observations. +// +// Adjust the interpolated background to satisfy mass continuity and +// match the supplied points exactly. +// +// Type: boolean + +adjust_background = FALSE; + +///////////// bkgd_obs_interpolation ////////////////// +// +// Interpolation method to fit background observations to the grid. +// +// TODO explain the various methods here. +// +// Type: enum +// Options: +// INTERP_NONE +// INTERP_SPLINE +// INTERP_KD_TREE +// INTERP_FRACTL + +bkgd_obs_interpolation = INTERP_SPLINE; + +//====================================================================== +// +// OPERATION SECTION. +// +//====================================================================== + +///////////// mode //////////////////////////////////// +// +// Coordinate system. +// +// XYZ for Cartesian, RTZ for spherical. +// +// Type: enum +// Options: +// MODE_XYZ +// MODE_RTZ + +mode = MODE_XYZ; + +///////////// mode //////////////////////////////////// +// Type of analysis. +// +// WIND for wind analysis +// THERMO for Thermo analysis using samurai output from a file +// WIND_THERMO for wind analysis and then Thermo analysis +// +// Type: enum +// Options: +// WIND +// THERMO +// WIND_THERMO + +analysis_type = THERMO; + +///////////// projection ////////////////////////////// +// +// Projection. +// +// TODO. +// +// Type: enum +// Options: +// PROJ_LAMBERT_CONFORMAL_CONIC +// PROJ_TRANSVERSE_MERCATOR_EXACT + +projection = PROJ_TRANSVERSE_MERCATOR_EXACT; + +///////////// data_directory ////////////////////////// +// +// Path to the data directory. +// +// Samurai will load data files from this directory. +// +// Type: string + +data_directory = "/glade/campaign/cisl/asap/samurai/data/beltrami"; + +///////////// output_directory //////////////////////// +// +// Path to the output directory. +// +// Samurai will write result files in this directory. +// +// Type: string + +output_directory = "."; + +///////////// preprocess_obs ////////////////////////// +// +// TODO. +// +// Type: boolean + +preprocess_obs = FALSE; + +///////////// num_iterations ////////////////////////// +// +// Max number of iterations to the multipass reduction factor. +// +// Multiple iterations will reduce the cutoff wavelengths and background +// error variance. +// +// Type: int + +num_iterations = 1; + +///////////// output_mish ///////////////////////////// +// +// Type: boolean + +output_mish = FALSE; + +///////////// output_txt ////////////////////////////// +// +// Type: boolean + +output_txt = FALSE; + +///////////// output_qc /////////////////////////////// +// +// Type: boolean + +output_qc = TRUE; + +///////////// output_netcdf /////////////////////////// +// +// Type: boolean + +output_netcdf = TRUE; + +///////////// output_asi ////////////////////////////// +// +// Type: boolean + +output_asi = FALSE; + +///////////// output_COAMPS /////////////////////////// +// +// Type: boolean + +output_COAMPS = FALSE; + +///////////// save_mish /////////////////////////////// +// +// Type: boolean + +save_mish = FALSE; + +//====================================================================== +// +// GRID DEFINITION SECTION. +// +//====================================================================== + +///////////// i_min /////////////////////////////////// +// +// Type: float + +i_min = 12.0; + +///////////// i_max /////////////////////////////////// +// +// Type: float + +i_max = 28.0; + +///////////// i_incr ////////////////////////////////// +// +// Type: float + +i_incr = 0.5; + +///////////// j_min /////////////////////////////////// +// +// Type: float + +j_min = -4.0; + +///////////// j_max /////////////////////////////////// +// +// Type: float + +j_max = 12.0; + +///////////// j_incr ////////////////////////////////// +// +// Type: float + +j_incr = 0.5; + +///////////// k_min /////////////////////////////////// +// +// Type: float + +k_min = 0.0; + +///////////// k_max /////////////////////////////////// +// +// Type: float + +k_max = 18.0; + +///////////// k_incr ////////////////////////////////// +// +// Type: float + +k_incr = 0.5; + +//====================================================================== +// +// BACKGROUND SECTION. +// +//====================================================================== + +///////////// ref_state /////////////////////////////// +// +// Type: string + +ref_state = "dunion_mt.snd"; + +///////////// ref_time //////////////////////////////// +// +// Reference time. +// +// hh:mm:ss. +// +// Type: string + +ref_time = "23:51:32"; + +///////////// i_background_roi //////////////////////// +// +// Radius of influence in the I direction. +// +// Type: float + +i_background_roi = 25; + +///////////// j_background_roi //////////////////////// +// +// Radius of influence in the J direction. +// +// Type: float + +j_background_roi = 25; + +//====================================================================== +// +// RADAR SECTION. +// +//====================================================================== + +///////////// radar_skip ////////////////////////////// +// +// Type: int + +radar_skip = 2; + +///////////// radar_stride //////////////////////////// +// +// Type: int + +radar_stride = 3; + +///////////// dynamic_stride ////////////////////////// +// +// Type: int + +dynamic_stride = 0; + +///////////// qr_variable ///////////////////////////// +// +// ???. +// +// Type: string + +qr_variable = "dbz"; + +///////////// radar_dbz /////////////////////////////// +// +// Radar reflectivity variable name in the input files. +// +// Example: DBZ for Eldora, DZ for CSU-CHILL. +// +// Type: string + +radar_dbz = "DBZ"; + +///////////// radar_vel /////////////////////////////// +// +// Radar velocity of scatterers away from instrument. +// +// Example: VG for Eldora, VE for CSU-CHILL, VR, ... +// +// Type: string + +radar_vel = "VG"; + +///////////// radar_sw //////////////////////////////// +// +// Type: string + +radar_sw = "SW"; + +///////////// i_reflectivity_roi ////////////////////// +// +// Type: float + +i_reflectivity_roi = 0.5; + +///////////// j_reflectivity_roi ////////////////////// +// +// Type: float + +j_reflectivity_roi = 0.5; + +///////////// k_reflectivity_roi ////////////////////// +// +// Type: float + +k_reflectivity_roi = 0.5; + +///////////// dbz_pseudow_weight ////////////////////// +// +// Type: float + +dbz_pseudow_weight = 1.0; + +///////////// mask_reflectivity /////////////////////// +// +// Type: float + +mask_reflectivity = -500; + +///////////// melting_zone_width ////////////////////// +// +// Type: float + +melting_zone_width = 1; + +///////////// mixed_phase_dbz ///////////////////////// +// +// Type: float + +mixed_phase_dbz = 20.0; + +///////////// rain_dbz //////////////////////////////// +// +// Type: float + +rain_dbz = 30.0; + +//====================================================================== +// +// PARAMETERS SECTION. +// +//====================================================================== + +//====================================================================== +// +// BOUNDARY CONDITIONS SECTION. +// +//====================================================================== + +///////////// i_rhou_bcL ////////////////////////////// +// +// Type: string + +i_rhou_bcL = "R0"; + +///////////// i_rhou_bcR ////////////////////////////// +// +// Type: string + +i_rhou_bcR = "R0"; + +///////////// i_rhov_bcL ////////////////////////////// +// +// Type: string + +i_rhov_bcL = "R0"; + +///////////// i_rhov_bcR ////////////////////////////// +// +// Type: string + +i_rhov_bcR = "R0"; + +///////////// i_rhow_bcL ////////////////////////////// +// +// Type: string + +i_rhow_bcL = "R0"; + +///////////// i_rhow_bcR ////////////////////////////// +// +// Type: string + +i_rhow_bcR = "R0"; + +///////////// i_tempk_bcL ///////////////////////////// +// +// Type: string + +i_tempk_bcL = "R0"; + +///////////// i_tempk_bcR ///////////////////////////// +// +// Type: string + +i_tempk_bcR = "R0"; + +///////////// i_qv_bcL //////////////////////////////// +// +// Type: string + +i_qv_bcL = "R0"; + +///////////// i_qv_bcR //////////////////////////////// +// +// Type: string + +i_qv_bcR = "R0"; + +///////////// i_rhoa_bcL ////////////////////////////// +// +// Type: string + +i_rhoa_bcL = "R0"; + +///////////// i_rhoa_bcR ////////////////////////////// +// +// Type: string + +i_rhoa_bcR = "R0"; + +///////////// i_qr_bcL //////////////////////////////// +// +// Type: string + +i_qr_bcL = "R0"; + +///////////// i_qr_bcR //////////////////////////////// +// +// Type: string + +i_qr_bcR = "R0"; + +///////////// j_rhou_bcL ////////////////////////////// +// +// Type: string + +j_rhou_bcL = "R0"; + +///////////// j_rhou_bcR ////////////////////////////// +// +// Type: string + +j_rhou_bcR = "R0"; + +///////////// j_rhov_bcL ////////////////////////////// +// +// Type: string + +j_rhov_bcL = "R0"; + +///////////// j_rhov_bcR ////////////////////////////// +// +// Type: string + +j_rhov_bcR = "R0"; + +///////////// j_rhow_bcL ////////////////////////////// +// +// Type: string + +j_rhow_bcL = "R0"; + +///////////// j_rhow_bcR ////////////////////////////// +// +// Type: string + +j_rhow_bcR = "R0"; + +///////////// j_tempk_bcL ///////////////////////////// +// +// Type: string + +j_tempk_bcL = "R0"; + +///////////// j_tempk_bcR ///////////////////////////// +// +// Type: string + +j_tempk_bcR = "R0"; + +///////////// j_qv_bcL //////////////////////////////// +// +// Type: string + +j_qv_bcL = "R0"; + +///////////// j_qv_bcR //////////////////////////////// +// +// Type: string + +j_qv_bcR = "R0"; + +///////////// j_rhoa_bcL ////////////////////////////// +// +// Type: string + +j_rhoa_bcL = "R0"; + +///////////// j_rhoa_bcR ////////////////////////////// +// +// Type: string + +j_rhoa_bcR = "R0"; + +///////////// j_qr_bcL //////////////////////////////// +// +// Type: string + +j_qr_bcL = "R0"; + +///////////// j_qr_bcR //////////////////////////////// +// +// Type: string + +j_qr_bcR = "R0"; + +///////////// k_rhou_bcL ////////////////////////////// +// +// Type: string + +k_rhou_bcL = "R0"; + +///////////// k_rhou_bcR ////////////////////////////// +// +// Type: string + +k_rhou_bcR = "R0"; + +///////////// k_rhov_bcL ////////////////////////////// +// +// Type: string + +k_rhov_bcL = "R0"; + +///////////// k_rhov_bcR ////////////////////////////// +// +// Type: string + +k_rhov_bcR = "R0"; + +///////////// k_rhow_bcL ////////////////////////////// +// +// Type: string + +k_rhow_bcL = "R1T0"; + +///////////// k_rhow_bcR ////////////////////////////// +// +// Type: string + +k_rhow_bcR = "R1T0"; + +///////////// k_tempk_bcL ///////////////////////////// +// +// Type: string + +k_tempk_bcL = "R0"; + +///////////// k_tempk_bcR ///////////////////////////// +// +// Type: string + +k_tempk_bcR = "R0"; + +///////////// k_qv_bcL //////////////////////////////// +// +// Type: string + +k_qv_bcL = "R0"; + +///////////// k_qv_bcR //////////////////////////////// +// +// Type: string + +k_qv_bcR = "R0"; + +///////////// k_rhoa_bcL ////////////////////////////// +// +// Type: string + +k_rhoa_bcL = "R0"; + +///////////// k_rhoa_bcR ////////////////////////////// +// +// Type: string + +k_rhoa_bcR = "R0"; + +///////////// k_qr_bcL //////////////////////////////// +// +// Type: string + +k_qr_bcL = "R0"; + +///////////// k_qr_bcR //////////////////////////////// +// +// Type: string + +k_qr_bcR = "R0"; + +//====================================================================== +// +// OBSERVATION ERRORS SECTION. +// +//====================================================================== + +///////////// dropsonde_rhoa_error //////////////////// +// +// Type: float + +dropsonde_rhoa_error = 1; + +///////////// dropsonde_rhou_error //////////////////// +// +// Type: float + +dropsonde_rhou_error = 2; + +///////////// dropsonde_rhov_error //////////////////// +// +// Type: float + +dropsonde_rhov_error = 2; + +///////////// dropsonde_rhow_error //////////////////// +// +// Type: float + +dropsonde_rhow_error = 50; + +///////////// dropsonde_tempk_error /////////////////// +// +// Type: float + +dropsonde_tempk_error = 2; + +///////////// dropsonde_qv_error ////////////////////// +// +// Type: float + +dropsonde_qv_error = 0.5; + +///////////// dropsonde_rhoua_error /////////////////// +// +// Type: float + +dropsonde_rhoua_error = 1; + +///////////// flightlevel_rhoa_error ////////////////// +// +// Type: float + +flightlevel_rhoa_error = 1; + +///////////// flightlevel_rhou_error ////////////////// +// +// Type: float + +flightlevel_rhou_error = 1; + +///////////// flightlevel_rhov_error ////////////////// +// +// Type: float + +flightlevel_rhov_error = 1; + +///////////// flightlevel_rhow_error ////////////////// +// +// Type: float + +flightlevel_rhow_error = 1; + +///////////// flightlevel_tempk_error ///////////////// +// +// Type: float + +flightlevel_tempk_error = 1; + +///////////// flightlevel_qv_error //////////////////// +// +// Type: float + +flightlevel_qv_error = 0.5; + +///////////// flightlevel_rhoua_error ///////////////// +// +// Type: float + +flightlevel_rhoua_error = 1; + +///////////// mtp_rhoa_error ////////////////////////// +// +// Type: float + +mtp_rhoa_error = 0; + +///////////// mtp_tempk_error ///////////////////////// +// +// Type: float + +mtp_tempk_error = 0; + +///////////// insitu_rhoa_error /////////////////////// +// +// Type: float + +insitu_rhoa_error = 1; + +///////////// insitu_rhou_error /////////////////////// +// +// Type: float + +insitu_rhou_error = 1; + +///////////// insitu_rhov_error /////////////////////// +// +// Type: float + +insitu_rhov_error = 1; + +///////////// insitu_rhow_error /////////////////////// +// +// Type: float + +insitu_rhow_error = 2; + +///////////// insitu_tempk_error ////////////////////// +// +// Type: float + +insitu_tempk_error = 1; + +///////////// insitu_qv_error ///////////////////////// +// +// Type: float + +insitu_qv_error = 0.5; + +///////////// insitu_rhoua_error ////////////////////// +// +// Type: float + +insitu_rhoua_error = 1; + +///////////// sfmr_windspeed_error //////////////////// +// +// Type: float + +sfmr_windspeed_error = 10; + +///////////// qscat_rhou_error //////////////////////// +// +// Type: float + +qscat_rhou_error = 2.5; + +///////////// qscat_rhov_error //////////////////////// +// +// Type: float + +qscat_rhov_error = 2.5; + +///////////// ascat_rhou_error //////////////////////// +// +// Type: float + +ascat_rhou_error = 2.5; + +///////////// ascat_rhov_error //////////////////////// +// +// Type: float + +ascat_rhov_error = 2.5; + +///////////// amv_rhou_error ////////////////////////// +// +// Type: float + +amv_rhou_error = 10; + +///////////// amv_rhov_error ////////////////////////// +// +// Type: float + +amv_rhov_error = 10; + +///////////// lidar_sw_error ////////////////////////// +// +// Type: float + +lidar_sw_error = 1; + +///////////// lidar_power_error /////////////////////// +// +// Type: float + +lidar_power_error = 50; + +///////////// lidar_min_error ///////////////////////// +// +// Type: float + +lidar_min_error = 1; + +///////////// radar_sw_error ////////////////////////// +// +// Type: float + +radar_sw_error = 1; + +///////////// radar_fallspeed_error /////////////////// +// +// Type: float + +radar_fallspeed_error = 2; + +///////////// radar_min_error ///////////////////////// +// +// Type: float + +radar_min_error = 1; + +///////////// aeri_qv_error /////////////////////////// +// +// Type: float + +aeri_qv_error = 1; + +///////////// aeri_rhoa_error ///////////////////////// +// +// Type: float + +aeri_rhoa_error = 1; + +///////////// aeri_rhou_error ///////////////////////// +// +// Type: float + +aeri_rhou_error = 1; + +///////////// aeri_rhov_error ///////////////////////// +// +// Type: float + +aeri_rhov_error = 1; + +///////////// aeri_rhow_error ///////////////////////// +// +// Type: float + +aeri_rhow_error = 1; + +///////////// aeri_tempk_error //////////////////////// +// +// Type: float + +aeri_tempk_error = 1; + +///////////// bg_obs_error //////////////////////////// +// +// Type: float + +bg_obs_error = 1; + +///////////// bg_interpolation_error ////////////////// +// +// Type: float + +bg_interpolation_error = 1; + +///////////// mesonet_qv_error //////////////////////// +// +// Type: float + +mesonet_qv_error = 1; + +///////////// mesonet_rhoa_error ////////////////////// +// +// Type: float + +mesonet_rhoa_error = 1; + +///////////// mesonet_rhou_error ////////////////////// +// +// Type: float + +mesonet_rhou_error = 1; + +///////////// mesonet_rhov_error ////////////////////// +// +// Type: float + +mesonet_rhov_error = 1; + +///////////// mesonet_rhow_error ////////////////////// +// +// Type: float + +mesonet_rhow_error = 1; + +///////////// mesonet_tempk_error ///////////////////// +// +// Type: float + +mesonet_tempk_error = 1; + +//====================================================================== +// +// XYP SPECIFIC SECTION. +// +//====================================================================== + +///////////// output_latlon_increment ///////////////// +// +// Type: float + +output_latlon_increment = -1; + +///////////// output_pressure_increment /////////////// +// +// Type: float + +output_pressure_increment = -1; + +//====================================================================== +// +// OPTION SECTION. +// +//====================================================================== + +///////////// max_radar_elevation ///////////////////// +// +// Type: float + +max_radar_elevation = 45; + +///////////// horizontal_radar_appx /////////////////// +// +// Type: boolean + +horizontal_radar_appx = TRUE; + +///////////// allow_background_missing_values ///////// +// +// Type: boolean + +allow_background_missing_values = TRUE; + +///////////// allow_negative_angles /////////////////// +// +// Type: boolean + +allow_negative_angles = TRUE; + +///////////// array_order ///////////////////////////// +// +// Order of compressed arrays passed through samurai.h interface. +// +// One of row_order or column_order. +// +// Type: string + +array_order = "row_order"; + +///////////// bg_interpolation //////////////////////// +// +// Type: string + +bg_interpolation = "none"; + +//====================================================================== +// +// KD TREE NEAREST NEIGHBOR SECTION. +// +//====================================================================== + +///////////// bkgd_kd_max_distance //////////////////// +// +// Any point outside of that distance will be ignored. +// +// Type: float + +bkgd_kd_max_distance = 20; + +///////////// bkgd_kd_num_neighbors /////////////////// +// +// Values will be averaged over these many nearest neighbors. +// +// Type: int + +bkgd_kd_num_neighbors = 1; + +//====================================================================== +// +// FRACTL INTERFACE SECTION. +// +//====================================================================== + +///////////// fractl_nc_file ////////////////////////// +// +// Need bkgd_obs_interpolation set to INTERP_FRACTL to be used. +// +// Type: string + +fractl_nc_file = ""; + +///////////// use_fractl_errors /////////////////////// +// +// Need bkgd_obs_interpolation set to 'fractl' to be used. +// +// Type: boolean + +use_fractl_errors = FALSE; + +//====================================================================== +// +// ITERATION DEPENDENT SECTION. +// +//====================================================================== + +///////////// mc_weight /////////////////////////////// +// +// Type: float +// 1D array - variable length. + +mc_weight = { 1, 1 }; + +///////////// neumann_u_weight //////////////////////// +// +// Type: float +// 1D array - variable length. + +neumann_u_weight = { 0.01, 0.01 }; + +///////////// neumann_v_weight //////////////////////// +// +// Type: float +// 1D array - variable length. + +neumann_v_weight = { 0.01, 0.01 }; + +///////////// dirichlet_w_weight ////////////////////// +// +// Type: float +// 1D array - variable length. + +dirichlet_w_weight = { 0.01, 0.01 }; + +///////////// bg_qr_error ///////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_qr_error = { 5, 1 }; + +///////////// bg_qv_error ///////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_qv_error = { 20, 2 }; + +///////////// bg_rhoa_error /////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_rhoa_error = { 20, 2 }; + +///////////// bg_rhou_error /////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_rhou_error = { 100, 2 }; + +///////////// bg_rhov_error /////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_rhov_error = { 100, 2 }; + +///////////// bg_rhow_error /////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_rhow_error = { 20, 2 }; + +///////////// bg_tempk_error ////////////////////////// +// +// Type: float +// 1D array - variable length. + +bg_tempk_error = { 20, 2 }; + +///////////// i_filter_length ///////////////////////// +// +// Type: float +// 1D array - variable length. + +i_filter_length = { 2, 2 }; + +///////////// j_filter_length ///////////////////////// +// +// Type: float +// 1D array - variable length. + +j_filter_length = { 2, 2 }; + +///////////// k_filter_length ///////////////////////// +// +// Type: float +// 1D array - variable length. + +k_filter_length = { 2, 2 }; + +///////////// i_spline_cutoff ///////////////////////// +// +// Type: float +// 1D array - variable length. + +i_spline_cutoff = { 2, 5 }; + +///////////// j_spline_cutoff ///////////////////////// +// +// Type: float +// 1D array - variable length. + +j_spline_cutoff = { 2, 5 }; + +///////////// k_spline_cutoff ///////////////////////// +// +// Type: float +// 1D array - variable length. + +k_spline_cutoff = { 2, 5 }; + +///////////// i_max_wavenumber //////////////////////// +// +// Type: float +// 1D array - variable length. + +i_max_wavenumber = { -1, -1 }; + +///////////// j_max_wavenumber //////////////////////// +// +// Type: float +// 1D array - variable length. + +j_max_wavenumber = { -1, -1 }; + +///////////// k_max_wavenumber //////////////////////// +// +// Type: float +// 1D array - variable length. + +k_max_wavenumber = { -1, -1 }; + diff --git a/ncar_scripts/TDRP/hurricane.tdrp b/ncar_scripts/TDRP/hurricane.tdrp index ed32e16..1316c32 100644 --- a/ncar_scripts/TDRP/hurricane.tdrp +++ b/ncar_scripts/TDRP/hurricane.tdrp @@ -147,6 +147,21 @@ bkgd_obs_interpolation = INTERP_SPLINE; mode = MODE_XYZ; +///////////// mode //////////////////////////////////// +// Type of analysis. +// +// WIND for wind analysis +// THERMO for Thermo analysis using samurai output from a file +// WIND_THERMO for wind analysis and then Thermo analysis +// +// Type: enum +// Options: +// WIND +// THERMO +// WIND_THERMO + +analysis_type = WIND; + ///////////// projection ////////////////////////////// // // Projection. diff --git a/ncar_scripts/TDRP/hurricane_4panel.tdrp b/ncar_scripts/TDRP/hurricane_4panel.tdrp index 30430a7..10fbd20 100644 --- a/ncar_scripts/TDRP/hurricane_4panel.tdrp +++ b/ncar_scripts/TDRP/hurricane_4panel.tdrp @@ -147,6 +147,21 @@ bkgd_obs_interpolation = INTERP_SPLINE; mode = MODE_XYZ; +///////////// mode //////////////////////////////////// +// Type of analysis. +// +// WIND for wind analysis +// THERMO for Thermo analysis using samurai output from a file +// WIND_THERMO for wind analysis and then Thermo analysis +// +// Type: enum +// Options: +// WIND +// THERMO +// WIND_THERMO + +analysis_type = WIND; + ///////////// projection ////////////////////////////// // // Projection. diff --git a/ncar_scripts/TDRP/supercell.tdrp b/ncar_scripts/TDRP/supercell.tdrp index 58bba8c..3940066 100644 --- a/ncar_scripts/TDRP/supercell.tdrp +++ b/ncar_scripts/TDRP/supercell.tdrp @@ -147,6 +147,21 @@ bkgd_obs_interpolation = INTERP_SPLINE; mode = MODE_XYZ; +///////////// mode //////////////////////////////////// +// Type of analysis. +// +// WIND for wind analysis +// THERMO for Thermo analysis using samurai output from a file +// WIND_THERMO for wind analysis and then Thermo analysis +// +// Type: enum +// Options: +// WIND +// THERMO +// WIND_THERMO + +analysis_type = WIND; + ///////////// projection ////////////////////////////// // // Projection. diff --git a/ncar_scripts/TDRP/typhoonChanthu2020.tdrp b/ncar_scripts/TDRP/typhoonChanthu2020.tdrp index f8168cd..62e94cc 100644 --- a/ncar_scripts/TDRP/typhoonChanthu2020.tdrp +++ b/ncar_scripts/TDRP/typhoonChanthu2020.tdrp @@ -169,6 +169,21 @@ bkgd_obs_interpolation = INTERP_SPLINE; mode = MODE_XYZ; +///////////// mode //////////////////////////////////// +// Type of analysis. +// +// WIND for wind analysis +// THERMO for Thermo analysis using samurai output from a file +// WIND_THERMO for wind analysis and then Thermo analysis +// +// Type: enum +// Options: +// WIND +// THERMO +// WIND_THERMO + +analysis_type = WIND; + ///////////// projection ////////////////////////////// // // Projection. diff --git a/src/Args.cpp b/src/Args.cpp index 82f48af..e2aafda 100644 --- a/src/Args.cpp +++ b/src/Args.cpp @@ -78,6 +78,10 @@ std::string interp_map[] = { "none", "spline", "kd_tree", "fractl" }; std::string mode_map[] = { "XYZ", "RTZ" }; +// This needs to match analysis_type_t in paramdef.samurai + +std::string analysis_type_map[] = { "WIND", "THERMO", "WIND_THERMO" }; + // This needs to match projection_t in paramdef.samurai std::string projection_map[] = { "lambert_conformal_conic", @@ -102,6 +106,7 @@ bool Args::paramsToHash(HashMap *configHash) { CONFIG_INSERT_BOOL(load_bg_coefficients); CONFIG_INSERT_FLOAT(mask_reflectivity); CONFIG_INSERT_MAP_VALUE(mode, mode_map); + CONFIG_INSERT_MAP_VALUE(analysis_type, analysis_type_map); CONFIG_INSERT_BOOL(output_asi); CONFIG_INSERT_BOOL(output_COAMPS); CONFIG_INSERT_STR(output_directory); diff --git a/src/Params.cc b/src/Params.cc index 219c4fa..5370d4f 100644 --- a/src/Params.cc +++ b/src/Params.cc @@ -741,6 +741,29 @@ 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->val_offset = (char *) &analysis_type - &_start_; + tt->enum_def.name = tdrpStrDup("analysis_type_t"); + tt->enum_def.nfields = 3; + tt->enum_def.fields = (enum_field_t *) + tdrpMalloc(tt->enum_def.nfields * sizeof(enum_field_t)); + tt->enum_def.fields[0].name = tdrpStrDup("WIND"); + tt->enum_def.fields[0].val = WIND; + tt->enum_def.fields[1].name = tdrpStrDup("THERMO"); + tt->enum_def.fields[1].val = THERMO; + tt->enum_def.fields[2].name = tdrpStrDup("WIND_THERMO"); + tt->enum_def.fields[2].val = WIND_THERMO; + tt->single_val.e = WIND; + tt++; + // Parameter 'projection' // ctype is '_projection_t' diff --git a/src/Params.hh b/src/Params.hh index bc334db..1fb2f93 100644 --- a/src/Params.hh +++ b/src/Params.hh @@ -79,6 +79,12 @@ public: MODE_RTZ = 1 } mode_t; + typedef enum { + WIND = 0, + THERMO = 1, + WIND_THERMO = 2 + } analysis_type_t; + typedef enum { PROJ_LAMBERT_CONFORMAL_CONIC = 0, PROJ_TRANSVERSE_MERCATOR_EXACT = 1 @@ -400,6 +406,8 @@ public: mode_t mode; + analysis_type_t analysis_type; + projection_t projection; char* data_directory; diff --git a/src/VarDriver3D.cpp b/src/VarDriver3D.cpp index 4dfdab1..4096417 100644 --- a/src/VarDriver3D.cpp +++ b/src/VarDriver3D.cpp @@ -101,6 +101,13 @@ bool VarDriver3D::validateDriver() 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") { + 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" ); diff --git a/src/paramdef.samurai b/src/paramdef.samurai index 9148e86..db401e2 100644 --- a/src/paramdef.samurai +++ b/src/paramdef.samurai @@ -104,6 +104,16 @@ paramdef enum mode_t { p_help = "XYZ for Cartesian, RTZ for spherical"; } mode; +typedef enum { + WIND, THERMO, WIND_THERMO +} analysis_type_t; + +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)" +} analysis_type; + typedef enum { PROJ_LAMBERT_CONFORMAL_CONIC, PROJ_TRANSVERSE_MERCATOR_EXACT diff --git a/src/thermo_reference/CostFunctionThermo.cpp b/src/thermo_reference/CostFunctionThermo.cpp new file mode 100644 index 0000000..4d72d1a --- /dev/null +++ b/src/thermo_reference/CostFunctionThermo.cpp @@ -0,0 +1,300 @@ +/* + * CostFunctionThermo.cpp + * samurai + * + * Copyright 2008 Michael Bell. All rights reserved. + * + */ + +#include "CostFunctionThermo.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +CostFunctionThermo::CostFunctionThermo(const int& numObs, const int& stateSize) + : CostFunction3D(numObs, stateSize) +{ +} + +CostFunctionThermo::~CostFunctionThermo() +{ +} + +void CostFunctionThermo::finalize() +{ +} + +void CostFunctionThermo::initialize(const QHash* config, real* bgU, real* obs, ReferenceState* ref) +{ + // Initialize number of variables + configHash = config; + + /* Set the output path */ + outputPath.setPath(configHash->value("output_directory")); + + // Horizontal boundary conditions + iBCL[0] = bcHash.value(configHash->value("i_pip_bcL")); + iBCR[0] = bcHash.value(configHash->value("i_pip_bcR")); + iBCL[1] = bcHash.value(configHash->value("i_thetarhop_bcL")); + iBCR[1] = bcHash.value(configHash->value("i_thetarhop_bcR")); + iBCL[2] = bcHash.value(configHash->value("i_ftheta_bcL")); + iBCR[2] = bcHash.value(configHash->value("i_ftheta_bcR")); + + jBCL[0] = bcHash.value(configHash->value("j_pip_bcL")); + jBCR[0] = bcHash.value(configHash->value("j_pip_bcR")); + jBCL[1] = bcHash.value(configHash->value("j_thetarhop_bcL")); + jBCR[1] = bcHash.value(configHash->value("j_thetarhop_bcR")); + jBCL[2] = bcHash.value(configHash->value("j_ftheta_bcL")); + jBCR[2] = bcHash.value(configHash->value("j_ftheta_bcR")); + + kBCL[0] = bcHash.value(configHash->value("k_pip_bcL")); + kBCR[0] = bcHash.value(configHash->value("k_pip_bcR")); + kBCL[1] = bcHash.value(configHash->value("k_thetarhop_bcL")); + kBCR[1] = bcHash.value(configHash->value("k_thetarhop_bcR")); + kBCL[2] = bcHash.value(configHash->value("k_ftheta_bcL")); + kBCR[2] = bcHash.value(configHash->value("k_ftheta_bcR")); + + // Define the Reference state + refstate = ref; + + // Assign local object pointers + bgFields = bgU; + rawObs = obs; + + iMin = configHash->value("i_min").toFloat(); + iMax = configHash->value("i_max").toFloat(); + DI = configHash->value("i_incr").toFloat(); + iDim = (int)((iMax - iMin)/DI) + 1; + jMin = configHash->value("j_min").toFloat(); + jMax = configHash->value("j_max").toFloat(); + DJ = configHash->value("j_incr").toFloat(); + jDim = (int)((jMax - jMin)/DJ) + 1; + kMin = configHash->value("k_min").toFloat(); + kMax = configHash->value("k_max").toFloat(); + DK = configHash->value("k_incr").toFloat(); + kDim = (int)((kMax - kMin)/DK) + 1; + + DIrecip = 1./DI; + DJrecip = 1./DJ; + DKrecip = 1./DK; + + // Adjust the internal, variable domain to include boundaries + adjustInternalDomain(1); + + // Define nodes with internal domain + int nodes = iDim*jDim*kDim; + + // Set up the initial recursive filter + iFilter = new RecursiveFilter(4); + jFilter = new RecursiveFilter(4); + kFilter = new RecursiveFilter(4); + + // Allocate memory for the needed arrays + // These are common to all CostFunctions + currState = new real[nState]; + currGradient = new real[nState]; + tempState = new real[nState]; + tempGradient = new real[nState]; + xt = new real[nState]; + df = new real[nState]; + + // These are local to this one + CTHTd = new real[nState]; + stateU = new real[nState]; + obsVector = new real[mObs*(obMetaSize+varDim*derivDim)]; + HCq = new real[mObs+nodes]; + innovation = new real[mObs+nodes]; + bgState = new real[nState]; + bgStdDev = new real[nState]; + stateA = new real[nState]; + stateB = new real[nState]; + stateC = new real[nState]; + + if (iBCL[0] == PERIODIC) { + iLDim = iDim-2; + } else { + iLDim = 4; + } + if (jBCL[0] == PERIODIC) { + jLDim = jDim-2; + } else { + jLDim = 4; + } + kLDim = 4; + for (int var = 0; var < varDim; ++var) { + iRank[var] = iDim - rankHash[iBCL[var]] - rankHash[iBCR[var]]; + if (iBCL[var] == PERIODIC) iRank[var]--; + jRank[var] = jDim - rankHash[jBCL[var]] - rankHash[jBCR[var]]; + if (jBCL[var] == PERIODIC) jRank[var]--; + kRank[var] = kDim - rankHash[kBCL[var]] - rankHash[kBCR[var]]; + // Need to verify that Rank is sufficient (>2?) + iL[var] = new real[iRank[var]*iLDim]; + jL[var] = new real[jRank[var]*jLDim]; + kL[var] = new real[kRank[var]*kLDim]; + iGamma[var] = new real[iRank[var]*iDim]; + jGamma[var] = new real[jRank[var]*jDim]; + kGamma[var] = new real[kRank[var]*kDim]; + + } + + // Precalculate the basis functions for lookup table option + /* + basisappx = configHash->value("spline_approximation").toInt(); + if (basisappx > 0) { + basis0 = new real[2000000]; + basis1 = new real[2000000]; + fillBasisLookup(); + } */ + + // Initialize the Fourier transforms + iFFTin = (double*) fftw_malloc(sizeof(double) * iDim); + iFFTout = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * iDim); + iForward = fftw_plan_dft_r2c_1d(iDim, iFFTin, iFFTout, FFTW_MEASURE); + iBackward = fftw_plan_dft_c2r_1d(iDim, iFFTout, iFFTin, FFTW_MEASURE); + jFFTin = (double*) fftw_malloc(sizeof(double) * jDim); + jFFTout = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * jDim); + jForward = fftw_plan_dft_r2c_1d(jDim, jFFTin, jFFTout, FFTW_MEASURE); + jBackward = fftw_plan_dft_c2r_1d(jDim, jFFTout, jFFTin, FFTW_MEASURE); + iMaxWavenumber = configHash->value("i_max_wavenumber").toFloat(); + jMaxWavenumber = configHash->value("j_max_wavenumber").toFloat(); +} + +void CostFunctionThermo::initState(const int iteration) +{ + // Clear the state vector + cout << "Initializing state vector..." << endl; + for (int n = 0; n < nState; n++) { + currState[n] = 0.0; + currGradient[n] = 0.0; + tempState[n] = 0.0; + tempGradient[n] = 0.0; + xt[n] = 0.0; + df[n] = 0.0; + stateA[n] = 0.0; + stateB[n] = 0.0; + stateC[n] = 0.0; + } + + + // Initialize background errors and filter scales + bgError[0] = configHash->value("bg_pip_error").toFloat(); + bgError[1] = configHash->value("bg_thetarhop_error").toFloat(); + bgError[2] = configHash->value("bg_ftheta_error").toFloat(); + bgError[3] = 0; + bgError[4] = 0; + bgError[5] = 0; + bgError[6] = 0; + + // Initialize filter scales + // Set up the recursive filter + iFilterScale = configHash->value("i_filter_length").toFloat(); + jFilterScale = configHash->value("j_filter_length").toFloat(); + kFilterScale = configHash->value("k_filter_length").toFloat(); + iFilter->setFilterLengthScale(iFilterScale); + jFilter->setFilterLengthScale(jFilterScale); + kFilter->setFilterLengthScale(kFilterScale); + + // Set up the spline matrices + setupSplines(); + + // Flag whether or not to print the subgrid information + if ((configHash->value("output_mish") == "true") or + (configHash->value("save_mish") == "true")) { + mishFlag = 1; + } else { + mishFlag = 0; + } + + // Mass continuity weight + mcWeight = configHash->value("mc_weight").toFloat(); + cout << "Mass continuity weight set to " << mcWeight << endl; + + if (iteration == 1) { + cout << "Initializing background..." << endl; + // Set up the background state + for (int n = 0; n < nState; n++) { + bgState[n] = 0.0; + // Initialize the std. dev. to 1 for the initial SC transform + bgStdDev[n] = 1.0; + } + + if (configHash->value("load_bg_coefficients") == "true") { + for (int n = 0; n < nState; n++) { + stateA[n] = bgFields[n]; + } + } else { + // SB Transform on the original bg fields + SBtransform(bgFields, stateB); + + // SA transform = bg B's -> bg A's + SAtransform(stateB, stateA); + } + + // FF transform to match background and increment + FFtransform(stateA, bgState); + + } + + for (int var = 0; var < varDim; var++) { + // Using a constant bg error variance for now, but this could be variable across the nodes + 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; + bgStdDev[bIndex] = bgError[var]; + } + } + } + } + + // Compute and display the variable BG errors and RMS of values + for (int var = 0; var < varDim; var++) { + 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; + 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; + } else { + cout << "Variable " << var << " RMS = " << varScale << "\t BG Error = " << bgError[var] + << " ( Infinite! %)" << endl; + } + } + + // Load the obs locally and weight the nonlinear observation operators by interpolated bg fields + obAdjustments(); + + // Calculate the H matrix operator + calcHmatrix(); + + // d = y - HXb + calcInnovation(); + + // Output the original background field + outputAnalysis("background", bgState); + + cout << "Beginning analysis...\n"; + + // HTd + calcHTranspose(innovation, stateC); + + FFtransform(stateC, stateA); + SAtransform(stateA, stateB); + SCtransform(stateB, CTHTd); + + //Htransform(stateB); +} diff --git a/src/thermo_reference/CostFunctionThermo.h b/src/thermo_reference/CostFunctionThermo.h new file mode 100644 index 0000000..54fdc1d --- /dev/null +++ b/src/thermo_reference/CostFunctionThermo.h @@ -0,0 +1,42 @@ +/* + * CostFunctionThermo.h + * samurai + * + * Copyright 2008 Michael Bell. All rights reserved. + * + */ + +#ifndef COSTFUNCTHERMO_H +#define COSTFUNCTHERMO_H + +#include "CostFunction3D.h" +#include "CostFunctionRTZ.h" +#include +#include +#include +#include +#include +#include +#include +#include +//#include +#include + +class CostFunctionThermo: public CostFunction3D +{ + +public: + CostFunctionThermo(const int& numObs = 0, const int& stateSize = 0); + ~CostFunctionThermo(); + void initialize(const QHash* config, real* bgU, real* obs, ReferenceState* ref); + void finalize(); + void initState(const int iteration); + +protected: + virtual bool outputAnalysis(const QString& suffix, real* Astate) = 0; + bool writeAsi(const QString& asiFileName); + bool writeNetCDF(const QString& netcdfFileName); + +}; + +#endif diff --git a/src/thermo_reference/CostFunctionThermoRTZ.cpp b/src/thermo_reference/CostFunctionThermoRTZ.cpp new file mode 100644 index 0000000..9fbb660 --- /dev/null +++ b/src/thermo_reference/CostFunctionThermoRTZ.cpp @@ -0,0 +1,573 @@ +#include "CostFunctionThermoRTZ.h" +#include +#include +#include +#include +#include +#include + +CostFunctionThermoRTZ::CostFunctionThermoRTZ(const int& numObs, const int& stateSize) + : CostFunctionThermo(numObs, stateSize) +{ +} + +CostFunctionThermoRTZ::~CostFunctionThermoRTZ() +{ +} + +bool CostFunctionThermoRTZ::outputAnalysis(const QString& suffix, real* Astate) +{ + + cout << "Outputting " << suffix.toStdString() << "...\n"; + // H --> to Mish for output + QString samuraiout = "samurai_RTZ_" + suffix + ".out"; + ofstream samuraistream; + if (configHash->value("output_txt") == "true") { + samuraistream.open(outputPath.absoluteFilePath(samuraiout).toAscii().data()); + samuraistream << "R\tT\tZ\tpip\tthetarhop\tftheta\n"; + samuraistream.precision(10); + } + + int analysisDim = 12; + int analysisSize = (iDim-2)*(jDim-2)*(kDim-2); + finalAnalysis = new real[analysisSize*analysisDim]; + + real Pi = acos(-1.0); + for (int iIndex = 1; iIndex < iDim-1; iIndex++) { + for (int ihalf = 0; ihalf <= mishFlag; ihalf++) { + for (int imu = -ihalf; imu <= ihalf; imu++) { + real i = iMin + DI * (iIndex + (0.5*sqrt(1./3.) * imu + 0.5*ihalf)); + real r = i*1000; + if (i > ((iDim-1)*DI + iMin)) continue; + + for (int jIndex = 1; jIndex < jDim-1; jIndex++) { + for (int jhalf =0; jhalf <=mishFlag; jhalf++) { + for (int jmu = -jhalf; jmu <= jhalf; jmu++) { + real j = jMin + DJ * (jIndex + (0.5*sqrt(1./3.) * jmu + 0.5*jhalf)); + if (j > ((jDim-1)*DJ + jMin)) continue; + + for (int kIndex = 1; kIndex < kDim-1; kIndex++) { + for (int khalf =0; khalf <=mishFlag; khalf++) { + for (int kmu = -khalf; kmu <= khalf; kmu++) { + real k = kMin + DK * (kIndex + (0.5*sqrt(1./3.) * kmu + 0.5*khalf)); + if (k > ((kDim-1)*DK + kMin)) continue; + + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0.; + real jbasis = 0.; + real kbasis = 0.; + real idbasis = 0.; + real jdbasis = 0.; + real kdbasis = 0.; + real pip = 0.; + real thetarhop = 0.; + real ftheta = 0.; + real pipdr = 0.; real thetarhopdr = 0.; real fthetadr = 0.; + real pipdt = 0.; real thetarhopdt = 0.; real fthetadt = 0.; + real pipdz = 0.; real thetarhopdz = 0.; real fthetadz = 0.; + + for (int var = 0; var < varDim; var++) { + for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 0, iBCL[var], iBCR[var]); + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 0, jBCL[var], jBCR[var]); + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 0, kBCL[var], kBCR[var]); + idbasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 1, iBCL[var], iBCR[var]); + jdbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 1, jBCL[var], jBCR[var]); + kdbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 1, kBCL[var], kBCR[var]); + real basis3x = ibasis*jbasis*kbasis; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + switch (var) { + case 0: + pip += Astate[aIndex] * basis3x; + pipdr += Astate[aIndex] * idbasis * jbasis * kbasis; + pipdt += 180 * Astate[aIndex] * ibasis * jdbasis * kbasis / (i * Pi); + pipdz += Astate[aIndex] * ibasis * jbasis * kdbasis; + break; + case 1: + thetarhop += Astate[aIndex + 1] * basis3x; + thetarhopdr += Astate[aIndex + 1] * idbasis * jbasis * kbasis; + thetarhopdt += 180 * Astate[aIndex + 1] * ibasis * jdbasis * kbasis / (i * Pi); + thetarhopdz += Astate[aIndex + 1] * ibasis * jbasis * kdbasis; + break; + case 2: + ftheta += Astate[aIndex + 2] * basis3x; + fthetadr += Astate[aIndex + 2] * idbasis * jbasis * kbasis; + fthetadt += 180 * Astate[aIndex + 2] * ibasis * jdbasis * kbasis / (i * Pi); + fthetadz += Astate[aIndex + 2] * ibasis * jbasis * kdbasis; + break; + } + } + } + } + } + + + // Avoid singularity at origin + if (r==0) { + pip = -999.; + thetarhop = -999.; + ftheta = -999.; + } + + + if (configHash->value("output_txt") == "true") { + samuraistream << scientific << i << "\t" << j << "\t" << k + << "\t" << pip << "\t" << thetarhop << "\t" << ftheta << "\n"; + } + + // On the nodes + if (!ihalf and !jhalf and !khalf){ + int fIndex = (iDim-2)*(jDim-2)*(kDim-2); + int posIndex = (iDim-2)*(jDim-2)*(kIndex-1) + (iDim-2)*(jIndex-1) + (iIndex-1); + finalAnalysis[fIndex * 0 + posIndex] = pip; + finalAnalysis[fIndex * 1 + posIndex] = thetarhop; + finalAnalysis[fIndex * 2 + posIndex] = ftheta; + finalAnalysis[fIndex * 3 + posIndex] = pipdr; + finalAnalysis[fIndex * 4 + posIndex] = thetarhopdr; + finalAnalysis[fIndex * 5 + posIndex] = fthetadr; + finalAnalysis[fIndex * 6 + posIndex] = pipdt; + finalAnalysis[fIndex * 7 + posIndex] = thetarhopdt; + finalAnalysis[fIndex * 8 + posIndex] = fthetadt; + finalAnalysis[fIndex * 9 + posIndex] = pipdz; + finalAnalysis[fIndex * 10 + posIndex] = thetarhopdz; + finalAnalysis[fIndex * 11 + posIndex] = fthetadz; + } + } + } + } + } + } + } + } + } + } + + QString fileName = "samurai_RTZ_" + suffix; + QString outFileName = outputPath.absoluteFilePath(fileName); + + // Write the Obs to a summary text file + if (configHash->value("output_qc") == "true") { + QString qcout = "samurai_QC_" + suffix + ".out"; + QString qcFileName = outputPath.absoluteFilePath(qcout); + ofstream qcstream(qcFileName.toAscii().data()); + ostream_iterator os(qcstream, "\t "); + *os++ = "Observation"; + *os++ = "Inverse Error"; + *os++ = "R"; + *os++ = "T"; + *os++ = "Z"; + *os++ = "Type"; + *os++ = "Time"; + *os++ = "pip"; + *os++ = "thetarhop"; + *os++ = "ftheta"; + *os++ = "Analysis"; + *os++ = "Background"; + qcstream << endl; + qcstream.precision(10); + + ostream_iterator od(qcstream, "\t "); + for (int m = 0; m < mObs; m++) { + int mi = m*(obMetaSize+varDim*derivDim); + real i = obsVector[mi+2]; + real j = obsVector[mi+3]; + real k = obsVector[mi+4]; + real tempsum = 0; + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0; + real jbasis = 0; + real kbasis = 0; + for (int var = 0; var < varDim; var++) { + for (int d = 0; d < derivDim; d++) { + int wgt_index = mi + obMetaSize +d*varDim + var; + if (!obsVector[wgt_index]) continue; + for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, derivative[d][2], kBCL[var], kBCR[var]); + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, derivative[d][1], iBCL[var], iBCR[var]); + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, derivative[d][0], jBCL[var], jBCR[var]); + tempsum += Astate[aIndex + var] * ibasis * jbasis * kbasis * obsVector[wgt_index]; + } + } + } + } + } + for (int t=0; t < obMetaSize-1; t++) { + *od++ = obsVector[mi+t]; + } + int unixtime = (int)obsVector[mi+6]; + QDateTime obtime; + obtime.setTime_t(unixtime); + obtime.setTimeSpec(Qt::UTC); + QString timestring = obtime.toString("hh:mm:ss.zzz"); + qcstream << timestring.toStdString() << "\t"; + + // Multiply the weight by the ob -- Observations.in has individual weights alreadt + // Only non-derivative for now + for (int t=obMetaSize; tvalue("output_netcdf") == "true") { + QString cdfFileName = outFileName + ".nc"; + if (!writeNetCDF(outputPath.absoluteFilePath(cdfFileName))) + cout << "Error writing netcdf file " << cdfFileName.toStdString() << endl; + } + // Write out to an asi file + if (configHash->value("output_asi") == "true") { + QString asiFileName = outFileName + ".asi"; + if (!writeAsi(outputPath.absoluteFilePath(asiFileName))) + cout << "Error writing asi file " << asiFileName.toStdString() << endl; + } + // Set the domain back + adjustInternalDomain(1); + + // Free the memory for the analysis variables + delete[] finalAnalysis; + + return true; + +} + +bool CostFunctionThermoRTZ::writeNetCDF(const QString& netcdfFileName) +{ + NcError err(NcError::verbose_nonfatal); + int NC_ERR = 0; + + // Create the file. + NcFile dataFile(netcdfFileName.toAscii(), NcFile::Replace); + + // Check to see if the file was created. + if(!dataFile.is_valid()) + return NC_ERR; + + // Define the dimensions. NetCDF will hand back an ncDim object for + // each. + NcDim *lvlDim, *radDim, *thetaDim, *timeDim; + if (!(radDim = dataFile.add_dim("radius", iDim))) + return NC_ERR; + if (!(thetaDim = dataFile.add_dim("theta", jDim))) + return NC_ERR; + if (!(lvlDim = dataFile.add_dim("altitude", kDim))) + return NC_ERR; + // Add an unlimited dimension... + if (!(timeDim = dataFile.add_dim("time"))) + return NC_ERR; + + // Define the coordinate variables. + NcVar *lvlVar, *timeVar, *radVar, *thetaVar; + /* NcVar *latVar, *lonVar; + if (!(lonVar = dataFile.add_var("longitude", ncFloat, radDim))) + return NC_ERR; + if (!(latVar = dataFile.add_var("latitude", ncFloat, thetaDim))) + return NC_ERR; */ + if (!(radVar = dataFile.add_var("radius", ncFloat, radDim))) + return NC_ERR; + if (!(thetaVar = dataFile.add_var("theta", ncFloat, thetaDim))) + return NC_ERR; + if (!(lvlVar = dataFile.add_var("altitude", ncFloat, lvlDim))) + return NC_ERR; + if (!(timeVar = dataFile.add_var("time", ncInt, timeDim))) + return NC_ERR; + + // Define units attributes for coordinate vars. This attaches a + // text attribute to each of the coordinate variables, containing + // the units. + /* if (!latVar->add_att("units", "degrees_north")) + return NC_ERR; + if (!lonVar->add_att("units", "degrees_east")) + return NC_ERR; */ + if (!radVar->add_att("units", "km")) + return NC_ERR; + if (!thetaVar->add_att("units", "degrees")) + return NC_ERR; + if (!lvlVar->add_att("units", "km")) + return NC_ERR; + if (!timeVar->add_att("units", "seconds since 1970-01-01 00:00:00 +0000")) + return NC_ERR; + + // Define the netCDF variables + NcVar *pip, *thetarhop, *ftheta; + NcVar *dpipdr, *dthetarhopdr, *dfthetadr, *dpipdt, *dthetarhopdt, *dfthetadt, *dpipdz, *dthetarhopdz, *dfthetadz; + + if (!(pip = dataFile.add_var("PIP", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(thetarhop = dataFile.add_var("THETARHOP", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(ftheta = dataFile.add_var("FTHETA", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + + return NC_ERR; + if (!(dpipdr = dataFile.add_var("DPIPDR", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dthetarhopdr = dataFile.add_var("DTHETARHOPDR", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dfthetadr = dataFile.add_var("DFTHETADR", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dpipdt = dataFile.add_var("DPIPDT", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dthetarhopdt = dataFile.add_var("DTHETARHOPDT", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dfthetadt = dataFile.add_var("DFTHETADT", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dpipdz = dataFile.add_var("DPIPDZ", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dthetarhopdz = dataFile.add_var("DTHETARHOPDZ", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + if (!(dfthetadz = dataFile.add_var("DFTHETADZ", ncFloat, timeDim, + lvlDim, thetaDim, radDim))) + return NC_ERR; + + + // Define units attributes for data variables. + if (!pip->add_att("units", "???")) + return NC_ERR; + if (!thetarhop->add_att("units", "???")) + return NC_ERR; + if (!ftheta->add_att("units", "???")) + return NC_ERR; + if (!dpipdr->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dthetarhopdr->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dfthetadr->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dpipdt->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dthetarhopdt->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dfthetadt->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dpipdz->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dthetarhopdz->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dfthetadz->add_att("units", "10-5s-1")) + return NC_ERR; + + // Define long names for data variables. + if (!pip->add_att("long_name", "pi prime")) + return NC_ERR; + if (!thetarhop->add_att("long_name", "theta rho prime")) + return NC_ERR; + if (!ftheta->add_att("long_name", "f theta")) + return NC_ERR; + if (!dpipdr->add_att("long_name", "pi prime gradient")) + return NC_ERR; + if (!dthetarhopdr->add_att("long_name", "theta rho prime gradient")) + return NC_ERR; + if (!dfthetadr->add_att("long_name", "f theta gradient")) + return NC_ERR; + if (!dpipdt->add_att("long_name", "pi prime gradient")) + return NC_ERR; + if (!dthetarhopdt->add_att("long_name", "theta rho prime gradient")) + return NC_ERR; + if (!dfthetadt->add_att("long_name", "f theta gradient")) + return NC_ERR; + if (!dpipdz->add_att("long_name", "pi prime gradient")) + return NC_ERR; + if (!dthetarhopdz->add_att("long_name", "theta rho prime gradient")) + return NC_ERR; + if (!dfthetadz->add_att("long_name", "f theta gradient")) + return NC_ERR; + + // Define missing data + if (!pip->add_att("missing_value", -999.f)) + return NC_ERR; + if (!thetarhop->add_att("missing_value", -999.f)) + return NC_ERR; + if (!ftheta->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dpipdr->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dthetarhopdr->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dfthetadr->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dpipdt->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dthetarhopdt->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dfthetadt->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dpipdz->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dthetarhopdz->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dfthetadz->add_att("missing_value", -999.f)) + return NC_ERR; + + // Define _Fill_Value for NCL users + if (!pip->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!thetarhop->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!ftheta->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dpipdr->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dthetarhopdr->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dfthetadr->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dpipdt->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dthetarhopdt->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dfthetadt->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dpipdz->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dthetarhopdz->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dfthetadz->add_att("_FillValue", -999.f)) + return NC_ERR; + + // Write the coordinate variable data to the file. + /* real *lons = new real[iDim]; + real *lats = new real[jDim]; */ + real *levs = new real[kDim]; + real *radius = new real[iDim]; + real *thetadeg = new real[jDim]; + int time[2]; + + // Reference time and position from center file + time[0] = configHash->value("ref_time").toInt(); + + /* real latReference = configHash->value("ref_lat").toFloat(); + real lonReference = configHash->value("ref_lon").toFloat(); + real refX, refY; + GeographicLib::TransverseMercatorExact tm = GeographicLib::TransverseMercatorExact::UTM(); + tm.Forward(lonReference, latReference, lonReference, refX, refY); */ + + for (int iIndex = 0; iIndex < iDim; iIndex++) { + real i = (iMin + DI * iIndex); + radius[iIndex] = i; + /* real j = (jMin + DJ * (jDim/2))*1000; + real latnull = 0; + tm.Reverse(lonReference,refX + i, refY + j, latnull, lons[iIndex]); + x[iIndex] = i/1000; */ + } + + for (int jIndex = 0; jIndex < jDim; jIndex++) { + real j = (jMin + DJ * jIndex); + thetadeg[jIndex] = j; + /* real i = (iMin + DI * (iDim/2))*1000; + real lonnull = 0; + tm.Reverse(lonReference,refX + i, refY + j, lats[jIndex], lonnull); + y[jIndex] = j/1000; */ + } + + /* if (!lonVar->put(lons, iDim)) + return NC_ERR; + + if (!latVar->put(lats, jDim)) + return NC_ERR; */ + + if (!radVar->put(radius, iDim)) + return NC_ERR; + + if (!thetaVar->put(thetadeg, jDim)) + return NC_ERR; + + for (int kIndex = 0; kIndex < kDim; kIndex++) { + real k = kMin + DK * kIndex; + levs[kIndex] = k; + } + if (!lvlVar->put(levs, kDim)) + return NC_ERR; + + if (!timeVar->put(time, 1)) + return NC_ERR; + + // Write the data. + for (int rec = 0; rec < 1; rec++) + { + if (!pip->put_rec(&finalAnalysis[0], rec)) + return NC_ERR; + if (!thetarhop->put_rec(&finalAnalysis[iDim*jDim*kDim*1], rec)) + return NC_ERR; + if (!ftheta->put_rec(&finalAnalysis[iDim*jDim*kDim*2], rec)) + return NC_ERR; + if (!dpipdr->put_rec(&finalAnalysis[iDim*jDim*kDim*3], rec)) + return NC_ERR; + if (!dthetarhopdr->put_rec(&finalAnalysis[iDim*jDim*kDim*4], rec)) + return NC_ERR; + if (!dfthetadr->put_rec(&finalAnalysis[iDim*jDim*kDim*5], rec)) + return NC_ERR; + if (!dpipdt->put_rec(&finalAnalysis[iDim*jDim*kDim*6], rec)) + return NC_ERR; + if (!dthetarhopdt->put_rec(&finalAnalysis[iDim*jDim*kDim*7], rec)) + return NC_ERR; + if (!dfthetadt->put_rec(&finalAnalysis[iDim*jDim*kDim*8], rec)) + return NC_ERR; + if (!dpipdz->put_rec(&finalAnalysis[iDim*jDim*kDim*9], rec)) + return NC_ERR; + if (!dthetarhopdz->put_rec(&finalAnalysis[iDim*jDim*kDim*10], rec)) + return NC_ERR; + if (!dfthetadz->put_rec(&finalAnalysis[iDim*jDim*kDim*11], rec)) + return NC_ERR; + } + + // The file is automatically closed by the destructor. This frees + // up any internal netCDF resources associated with the file, and + // flushes any buffers. + /* delete[] lats; + delete[] lons; */ + delete[] levs; + delete[] radius; + delete[] thetadeg; + return true; + +} + +bool CostFunctionThermoRTZ::writeAsi(const QString& asiFileName) +{ + cout << "Writing Asi not implemented yet" << endl; + + return true; +} diff --git a/src/thermo_reference/CostFunctionThermoRTZ.h b/src/thermo_reference/CostFunctionThermoRTZ.h new file mode 100644 index 0000000..c147b41 --- /dev/null +++ b/src/thermo_reference/CostFunctionThermoRTZ.h @@ -0,0 +1,19 @@ +#ifndef COSTFUNCTIONTHERMORTZ_H +#define COSTFUNCTIONTHERMORTZ_H + +#include "CostFunctionThermo.h" +#include + +class CostFunctionThermoRTZ: public CostFunctionThermo +{ +public: + CostFunctionThermoRTZ(const int& numObs = 0, const int& stateSize = 0); + ~CostFunctionThermoRTZ(); + +private: + bool outputAnalysis(const QString& suffix, real* Astate); + bool writeAsi(const QString& asiFileName); + bool writeNetCDF(const QString& netcdfFileName); +}; + +#endif // COSTFUNCTIONTHERMORTZ_H diff --git a/src/thermo_reference/CostFunctionThermoXYZ.cpp b/src/thermo_reference/CostFunctionThermoXYZ.cpp new file mode 100644 index 0000000..9cb0d5e --- /dev/null +++ b/src/thermo_reference/CostFunctionThermoXYZ.cpp @@ -0,0 +1,559 @@ +#include "CostFunctionThermoXYZ.h" +#include +#include +#include +#include +#include +#include + +CostFunctionThermoXYZ::CostFunctionThermoXYZ(const int& numObs, const int& stateSize) + : CostFunctionThermo(numObs, stateSize) +{ +} + +CostFunctionThermoXYZ::~CostFunctionThermoXYZ() +{ +} + +bool CostFunctionThermoXYZ::outputAnalysis(const QString& suffix, real* Astate) +{ + + cout << "Outputting " << suffix.toStdString() << "...\n"; + // H --> to Mish for output + QString samuraiout = "samurai_XYZ_" + suffix + ".out"; + ofstream samuraistream; + if (configHash->value("output_txt") == "true") { + samuraistream.open(outputPath.absoluteFilePath(samuraiout).toAscii().data()); + samuraistream << "X\tY\tZ\tpip\tthetarhop\tftheta\n"; + samuraistream.precision(10); + } + + int analysisDim = 12; + int analysisSize = (iDim-2)*(jDim-2)*(kDim-2); + finalAnalysis = new real[analysisSize*analysisDim]; + real gausspoint = 0.5*sqrt(1./3.); + + for (int iIndex = 1; iIndex < iDim-1; iIndex++) { + for (int ihalf = 0; ihalf <= mishFlag; ihalf++) { + for (int imu = -ihalf; imu <= ihalf; imu++) { + real i = iMin + DI * (iIndex + (gausspoint * imu + 0.5*ihalf)); + if (i > ((iDim-1)*DI + iMin)) continue; + + for (int jIndex = 1; jIndex < jDim-1; jIndex++) { + for (int jhalf =0; jhalf <=mishFlag; jhalf++) { + for (int jmu = -jhalf; jmu <= jhalf; jmu++) { + real j = jMin + DJ * (jIndex + (gausspoint * jmu + 0.5*jhalf)); + if (j > ((jDim-1)*DJ + jMin)) continue; + + for (int kIndex = 1; kIndex < kDim-1; kIndex++) { + for (int khalf =0; khalf <=mishFlag; khalf++) { + for (int kmu = -khalf; kmu <= khalf; kmu++) { + real k = kMin + DK * (kIndex + (gausspoint * kmu + 0.5*khalf)); + if (k > ((kDim-1)*DK + kMin)) continue; + + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0.; + real jbasis = 0.; + real kbasis = 0.; + real idbasis = 0.; + real jdbasis = 0.; + real kdbasis = 0.; + real pip = 0.; + real thetarhop = 0.; + real ftheta = 0.; + real pipdx = 0.; real thetarhopdx = 0.; real fthetadx = 0.; + real pipdy = 0.; real thetarhopdy = 0.; real fthetady = 0.; + real pipdz = 0.; real thetarhopdz = 0.; real fthetadz = 0.; + + for (int var = 0; var < varDim; var++) { + for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 0, iBCL[var], iBCR[var]); + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 0, jBCL[var], jBCR[var]); + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 0, kBCL[var], kBCR[var]); + idbasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, 1, iBCL[var], iBCR[var]); + jdbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, 1, jBCL[var], jBCR[var]); + kdbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, 1, kBCL[var], kBCR[var]); + real basis3x = ibasis*jbasis*kbasis; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + switch (var) { + case 0: + pip += Astate[aIndex] * basis3x; + pipdx += Astate[aIndex] * idbasis * jbasis * kbasis; + pipdy += Astate[aIndex] * ibasis * jdbasis * kbasis; + pipdz += Astate[aIndex] * ibasis * jbasis * kdbasis; + break; + case 1: + thetarhop += Astate[aIndex + 1] * basis3x; + thetarhopdx += Astate[aIndex + 1] * idbasis * jbasis * kbasis; + thetarhopdy += Astate[aIndex + 1] * ibasis * jdbasis * kbasis; + thetarhopdz += Astate[aIndex + 1] * ibasis * jbasis * kdbasis; + break; + case 2: + ftheta += Astate[aIndex + 2] * basis3x; + fthetadx += Astate[aIndex + 2] * idbasis * jbasis * kbasis; + fthetady += Astate[aIndex + 2] * ibasis * jdbasis * kbasis; + fthetadz += Astate[aIndex + 2] * ibasis * jbasis * kdbasis; + break; + } + } + } + } + } + + + if (configHash->value("output_txt") == "true") { + samuraistream << scientific << i << "\t" << j << "\t" << k + << "\t" << pip << "\t" << thetarhop << "\t" << ftheta << "\t" << "\n"; + } + + // On the nodes + if (!ihalf and !jhalf and !khalf){ + int fIndex = (iDim-2)*(jDim-2)*(kDim-2); + int posIndex = (iDim-2)*(jDim-2)*(kIndex-1) + (iDim-2)*(jIndex-1) + (iIndex-1); + finalAnalysis[fIndex * 0 + posIndex] = pip; + finalAnalysis[fIndex * 1 + posIndex] = thetarhop; + finalAnalysis[fIndex * 2 + posIndex] = ftheta; + finalAnalysis[fIndex * 3 + posIndex] = pipdx; + finalAnalysis[fIndex * 4 + posIndex] = thetarhopdx; + finalAnalysis[fIndex * 5 + posIndex] = fthetadx; + finalAnalysis[fIndex * 6 + posIndex] = pipdy; + finalAnalysis[fIndex * 7 + posIndex] = thetarhopdy; + finalAnalysis[fIndex * 8 + posIndex] = fthetady; + finalAnalysis[fIndex * 9 + posIndex] = pipdz; + finalAnalysis[fIndex * 10 + posIndex] = thetarhopdz; + finalAnalysis[fIndex * 11 + posIndex] = fthetadz; + } + } + } + } + } + } + } + } + } + } + + QString fileName = "samurai_XYZ_" + suffix; + QString outFileName = outputPath.absoluteFilePath(fileName); + + // Write the Obs to a summary text file + if (configHash->value("output_qc") == "true") { + QString qcout = "samurai_QC_" + suffix + ".out"; + QString qcFileName = outputPath.absoluteFilePath(qcout); + ofstream qcstream(qcFileName.toAscii().data()); + ostream_iterator os(qcstream, "\t "); + *os++ = "Observation"; + *os++ = "Inverse Error"; + *os++ = "X"; + *os++ = "Y"; + *os++ = "Z"; + *os++ = "Type"; + *os++ = "Time"; + *os++ = "pip"; + *os++ = "thetarhop"; + *os++ = "ftheta"; + *os++ = "Analysis"; + *os++ = "Background"; + qcstream << endl; + qcstream.precision(10); + + ostream_iterator od(qcstream, "\t "); + for (int m = 0; m < mObs; m++) { + int mi = m*(obMetaSize+varDim*derivDim); + real i = obsVector[mi+2]; + real j = obsVector[mi+3]; + real k = obsVector[mi+4]; + real tempsum = 0; + int ii = (int)((i - iMin)*DIrecip); + int jj = (int)((j - jMin)*DJrecip); + int kk = (int)((k - kMin)*DKrecip); + real ibasis = 0; + real jbasis = 0; + real kbasis = 0; + for (int var = 0; var < varDim; var++) { + for (int d = 0; d < derivDim; d++) { + int wgt_index = mi + obMetaSize +d*varDim + var; + if (!obsVector[wgt_index]) continue; + for (int kNode = max(kk-1,0); kNode <= min(kk+2,kDim-1); ++kNode) { + kbasis = Basis(kNode, k, kDim-1, kMin, DK, DKrecip, derivative[d][2], kBCL[var], kBCR[var]); + for (int iiNode = (ii-1); iiNode <= (ii+2); ++iiNode) { + int iNode = iiNode; + if ((iBCL[var] == PERIODIC) and (iNode < 1)) iNode = iDim-3; + if ((iBCR[var] == PERIODIC) and (iNode > (iDim-3))) iNode = iiNode - (iDim-3); + if ((iNode < 0) or (iNode >= iDim)) continue; + ibasis = Basis(iNode, i, iDim-1, iMin, DI, DIrecip, derivative[d][1], iBCL[var], iBCR[var]); + + for (int jjNode = (jj-1); jjNode <= (jj+2); ++jjNode) { + int jNode = jjNode; + if ((jBCL[var] == PERIODIC) and (jNode < 1)) jNode = jDim-3; + if ((jBCR[var] == PERIODIC) and (jNode > (jDim-3))) jNode = jjNode - (jDim-3); + if ((jNode < 0) or (jNode >= jDim)) continue; + int aIndex = varDim*iDim*jDim*kNode + varDim*iDim*jNode +varDim*iNode; + jbasis = Basis(jNode, j, jDim-1, jMin, DJ, DJrecip, derivative[d][0], jBCL[var], jBCR[var]); + tempsum += Astate[aIndex + var] * ibasis * jbasis * kbasis * obsVector[wgt_index]; + } + } + } + } + } + for (int t=0; t < obMetaSize-1; t++) { + *od++ = obsVector[mi+t]; + } + int unixtime = (int)obsVector[mi+6]; + QDateTime obtime; + obtime.setTime_t(unixtime); + obtime.setTimeSpec(Qt::UTC); + QString timestring = obtime.toString("hh:mm:ss.zzz"); + qcstream << timestring.toStdString() << "\t"; + + // Multiply the weight by the ob -- Observations.in has individual weights already + // Only non-derivative for now + for (int t=obMetaSize; tvalue("output_netcdf") == "true") { + QString cdfFileName = outFileName + ".nc"; + if (!writeNetCDF(outputPath.absoluteFilePath(cdfFileName))) + cout << "Error writing netcdf file " << cdfFileName.toStdString() << endl; + } + // Write out to an asi file + if (configHash->value("output_asi") == "true") { + QString asiFileName = outFileName + ".asi"; + if (!writeAsi(outputPath.absoluteFilePath(asiFileName))) + cout << "Error writing asi file " << asiFileName.toStdString() << endl; + } + // Set the domain back + adjustInternalDomain(1); + + // Free the memory for the analysis variables + delete[] finalAnalysis; + + return true; + +} + +bool CostFunctionThermoXYZ::writeNetCDF(const QString& netcdfFileName) +{ + NcError err(NcError::verbose_nonfatal); + int NC_ERR = 0; + + // Create the file. + NcFile dataFile(netcdfFileName.toAscii(), NcFile::Replace); + + // Check to see if the file was created. + if(!dataFile.is_valid()) + return NC_ERR; + + // Define the dimensions. NetCDF will hand back an ncDim object for + // each. + NcDim *lvlDim, *latDim, *lonDim, *timeDim; + if (!(lonDim = dataFile.add_dim("longitude", iDim))) + return NC_ERR; + if (!(latDim = dataFile.add_dim("latitude", jDim))) + return NC_ERR; + if (!(lvlDim = dataFile.add_dim("altitude", kDim))) + return NC_ERR; + // Add an unlimited dimension... + if (!(timeDim = dataFile.add_dim("time"))) + return NC_ERR; + + // Define the coordinate variables. + NcVar *latVar, *lonVar, *lvlVar, *timeVar, *xVar, *yVar; + if (!(lonVar = dataFile.add_var("longitude", ncFloat, lonDim))) + return NC_ERR; + if (!(latVar = dataFile.add_var("latitude", ncFloat, latDim))) + return NC_ERR; + if (!(xVar = dataFile.add_var("x", ncFloat, lonDim))) + return NC_ERR; + if (!(yVar = dataFile.add_var("y", ncFloat, latDim))) + return NC_ERR; + if (!(lvlVar = dataFile.add_var("altitude", ncFloat, lvlDim))) + return NC_ERR; + if (!(timeVar = dataFile.add_var("time", ncInt, timeDim))) + return NC_ERR; + + // Define units attributes for coordinate vars. This attaches a + // text attribute to each of the coordinate variables, containing + // the units. + if (!latVar->add_att("units", "degrees_north")) + return NC_ERR; + if (!lonVar->add_att("units", "degrees_east")) + return NC_ERR; + if (!xVar->add_att("units", "km")) + return NC_ERR; + if (!yVar->add_att("units", "km")) + return NC_ERR; + if (!lvlVar->add_att("units", "km")) + return NC_ERR; + if (!timeVar->add_att("units", "seconds since 1970-01-01 00:00:00 +0000")) + return NC_ERR; + + // Define the netCDF variables + NcVar *pip, *thetarhop, *ftheta; + NcVar *dpipdx, *dthetarhopdx, *dfthetadx, *dpipdy, *dthetarhopdy, *dfthetady, *dpipdz, *dthetarhopdz, *dfthetadz; + + + if (!(pip = dataFile.add_var("PIP", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(thetarhop = dataFile.add_var("THETARHOP", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(ftheta = dataFile.add_var("FTHETA", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dpipdx = dataFile.add_var("DPIPDX", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dthetarhopdx = dataFile.add_var("DTHETARHOPDX", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dfthetadx = dataFile.add_var("DFTHETADX", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dpipdy = dataFile.add_var("DPIPDY", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dthetarhopdy = dataFile.add_var("DTHETARHOPDY", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dfthetady = dataFile.add_var("DFTHETADY", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dpipdz = dataFile.add_var("DPIPDZ", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dthetarhopdz = dataFile.add_var("DTHETARHOPDZ", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + if (!(dfthetadz = dataFile.add_var("DFTHETADZ", ncFloat, timeDim, + lvlDim, latDim, lonDim))) + return NC_ERR; + + // Define units attributes for data variables. + if (!pip->add_att("units", "???")) + return NC_ERR; + if (!thetarhop->add_att("units", "???")) + return NC_ERR; + if (!ftheta->add_att("units", "???")) + return NC_ERR; + if (!dpipdx->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dthetarhopdx->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dfthetadx->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dpipdy->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dthetarhopdy->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dfthetady->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dpipdz->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dthetarhopdz->add_att("units", "10-5s-1")) + return NC_ERR; + if (!dfthetadz->add_att("units", "10-5s-1")) + return NC_ERR; + + // Define long names for data variables. + if (!pip->add_att("long_name", "pi prime")) + return NC_ERR; + if (!thetarhop->add_att("long_name", "theta rho prime")) + return NC_ERR; + if (!ftheta->add_att("long_name", "f theta")) + return NC_ERR; + if (!dpipdx->add_att("long_name", "pi prime gradient")) + return NC_ERR; + if (!dthetarhopdx->add_att("long_name", "theta rho prime gradient")) + return NC_ERR; + if (!dfthetadx->add_att("long_name", "f theta gradient")) + return NC_ERR; + if (!dpipdy->add_att("long_name", "pi prime gradient")) + return NC_ERR; + if (!dthetarhopdy->add_att("long_name", "theta rho prime gradient")) + return NC_ERR; + if (!dfthetady->add_att("long_name", "f theta gradient")) + return NC_ERR; + if (!dpipdz->add_att("long_name", "pi prime gradient")) + return NC_ERR; + if (!dthetarhopdz->add_att("long_name", "theta rho prime gradient")) + return NC_ERR; + if (!dfthetadz->add_att("long_name", "f theta gradient")) + return NC_ERR; + + // Define missing data + if (!pip->add_att("missing_value", -999.f)) + return NC_ERR; + if (!thetarhop->add_att("missing_value", -999.f)) + return NC_ERR; + if (!ftheta->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dpipdx->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dthetarhopdx->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dfthetadx->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dpipdy->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dthetarhopdy->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dfthetady->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dpipdz->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dthetarhopdz->add_att("missing_value", -999.f)) + return NC_ERR; + if (!dfthetadz->add_att("missing_value", -999.f)) + return NC_ERR; + + // Define _Fill_Value for NCL users + if (!pip->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!thetarhop->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!ftheta->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dpipdx->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dthetarhopdx->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dfthetadx->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dpipdy->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dthetarhopdy->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dfthetady->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dpipdz->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dthetarhopdz->add_att("_FillValue", -999.f)) + return NC_ERR; + if (!dfthetadz->add_att("_FillValue", -999.f)) + return NC_ERR; + + // Write the coordinate variable data to the file. + real *lons = new real[iDim]; + real *lats = new real[jDim]; + real *levs = new real[kDim]; + real *x = new real[iDim]; + real *y = new real[jDim]; + int time[2]; + + // Reference time and position from center file + time[0] = configHash->value("ref_time").toInt(); + real latReference = configHash->value("ref_lat").toFloat(); + real lonReference = configHash->value("ref_lon").toFloat(); + real refX, refY; + + GeographicLib::TransverseMercatorExact tm = GeographicLib::TransverseMercatorExact::UTM(); + tm.Forward(lonReference, latReference, lonReference, refX, refY); + for (int iIndex = 0; iIndex < iDim; iIndex++) { + real i = (iMin + DI * iIndex)*1000; + real j = (jMin + DJ * (jDim/2))*1000; + real latnull = 0; + tm.Reverse(lonReference,refX + i, refY + j, latnull, lons[iIndex]); + x[iIndex] = i/1000; + } + + for (int jIndex = 0; jIndex < jDim; jIndex++) { + real i = (iMin + DI * (iDim/2))*1000; + real j = (jMin + DJ * jIndex)*1000; + real lonnull = 0; + tm.Reverse(lonReference,refX + i, refY + j, lats[jIndex], lonnull); + y[jIndex] = j/1000; + } + + if (!lonVar->put(lons, iDim)) + return NC_ERR; + + if (!latVar->put(lats, jDim)) + return NC_ERR; + + if (!xVar->put(x, iDim)) + return NC_ERR; + + if (!yVar->put(y, jDim)) + return NC_ERR; + + for (int kIndex = 0; kIndex < kDim; kIndex++) { + real k = kMin + DK * kIndex; + levs[kIndex] = k; + } + if (!lvlVar->put(levs, kDim)) + return NC_ERR; + + if (!timeVar->put(time, 1)) + return NC_ERR; + + // Write the data. + for (int rec = 0; rec < 1; rec++) + { + if (!pip->put_rec(&finalAnalysis[0], rec)) + return NC_ERR; + if (!thetarhop->put_rec(&finalAnalysis[iDim*jDim*kDim*1], rec)) + return NC_ERR; + if (!ftheta->put_rec(&finalAnalysis[iDim*jDim*kDim*2], rec)) + return NC_ERR; + if (!dpipdx->put_rec(&finalAnalysis[iDim*jDim*kDim*3], rec)) + return NC_ERR; + if (!dthetarhopdx->put_rec(&finalAnalysis[iDim*jDim*kDim*4], rec)) + return NC_ERR; + if (!dfthetadx->put_rec(&finalAnalysis[iDim*jDim*kDim*5], rec)) + return NC_ERR; + if (!dpipdy->put_rec(&finalAnalysis[iDim*jDim*kDim*6], rec)) + return NC_ERR; + if (!dthetarhopdy->put_rec(&finalAnalysis[iDim*jDim*kDim*7], rec)) + return NC_ERR; + if (!dfthetady->put_rec(&finalAnalysis[iDim*jDim*kDim*8], rec)) + return NC_ERR; + if (!dpipdz->put_rec(&finalAnalysis[iDim*jDim*kDim*9], rec)) + return NC_ERR; + if (!dthetarhopdz->put_rec(&finalAnalysis[iDim*jDim*kDim*10], rec)) + return NC_ERR; + if (!dfthetadz->put_rec(&finalAnalysis[iDim*jDim*kDim*11], rec)) + return NC_ERR; + } + + // The file is automatically closed by the destructor. This frees + // up any internal netCDF resources associated with the file, and + // flushes any buffers. + delete[] lats; + delete[] lons; + delete[] levs; + delete[] x; + delete[] y; + return true; + +} + +bool CostFunctionThermoXYZ::writeAsi(const QString& asiFileName) +{ + cout << "Writing Asi not implemented yet" << endl; + return true; +} diff --git a/src/thermo_reference/CostFunctionThermoXYZ.h b/src/thermo_reference/CostFunctionThermoXYZ.h new file mode 100644 index 0000000..479dd1c --- /dev/null +++ b/src/thermo_reference/CostFunctionThermoXYZ.h @@ -0,0 +1,24 @@ +#ifndef COSTFUNCTIONTHERMOXYZ_H +#define COSTFUNCTIONTHERMOXYZ_H + +#include "CostFunctionThermo.h" +#include + +class CostFunctionThermoXYZ: public CostFunctionThermo +{ +public: + CostFunctionThermoXYZ(const int& numObs = 0, const int& stateSize = 0); + ~CostFunctionThermoXYZ(); + +private: + bool outputAnalysis(const QString& suffix, real* Astate); + bool writeAsi(const QString& asiFileName); + bool writeNetCDF(const QString& netcdfFileName); +}; + +#endif // COSTFUNCTIONTHERMOXYZ_H + + + + + diff --git a/src/thermo_reference/NetCDF.cpp b/src/thermo_reference/NetCDF.cpp new file mode 100644 index 0000000..baf6b6c --- /dev/null +++ b/src/thermo_reference/NetCDF.cpp @@ -0,0 +1,34 @@ +/* + * NetCDF.cpp + * samurai + * + * Created by Annette Foerster on 8/15/13. + * Adapted from $Id: pres_temp_4D_rd.cpp,v 1.11 2006/08/22 19:22:06 ed Exp $ + * (http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial) + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + + +#include "NetCDF.h" +#include +#include + +NetCDF::NetCDF(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim) :c_p(1005.7), g(9.81), f(5.85e-05), pi(3.14159265358979323846) + { + NDIMS = 4; + NALT = metFile_kdim; + NX = metFile_idim; // NRADIUS or NLON + NY = metFile_jdim; // NTHETA or NLAT + NREC = 0; + NC_ERR = 2; + + std::cout << "Constructor \n"; +} + +NetCDF::~NetCDF() { +} + + + + diff --git a/src/thermo_reference/NetCDF.h b/src/thermo_reference/NetCDF.h new file mode 100644 index 0000000..3631f6d --- /dev/null +++ b/src/thermo_reference/NetCDF.h @@ -0,0 +1,44 @@ + /* + * NetCDF.h + * samurai + * + * Created by Annette Foerster on 8/15/13. + * Based on example code from netCDF library + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + +#ifndef NETCDF_H +#define NETCDF_H + +#include +#include + +class NetCDF +{ + +public: + NetCDF(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim); + virtual ~NetCDF(); + + virtual int readNetCDF(const char* filename)=0; + virtual double getValue(const int &i,const int &j,const int &k,const QString &varName)=0; + virtual double getDerivative(const int &i,const int &j,const int &k, const QString &var, const int &der)=0; + virtual double calc_A(const int &i,const int &j,const int &k)=0; + virtual double calc_B(const int &i,const int &j,const int &k)=0; + virtual double calc_C(const int &i,const int &j,const int &k)=0; + virtual double calc_D(const int &i,const int &j,const int &k)=0; + virtual double calc_E(const int &i,const int &j,const int &k)=0; + + +protected: + int NDIMS, NALT, NX, NY, NREC, NC_ERR; + + const float c_p; + const float g; + const double f; //this needs to be made dynamic eventually, here lat assumed to be 22 deg north + const double pi; +}; + + +#endif diff --git a/src/thermo_reference/NetCDF_XYZ.cpp b/src/thermo_reference/NetCDF_XYZ.cpp new file mode 100644 index 0000000..c59cae1 --- /dev/null +++ b/src/thermo_reference/NetCDF_XYZ.cpp @@ -0,0 +1,449 @@ +/* + * NetCDF_RTZ.cpp + * samurai + * + * Created by Annette Foerster. + * Adapted from $Id: pres_temp_4D_rd.cpp,v 1.11 2006/08/22 19:22:06 ed Exp $ + * (http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial) + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + +#include "NetCDF.h" +#include "NetCDF_XYZ.h" +#include +#include +#include + + +NetCDF_XYZ::NetCDF_XYZ(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim) +:NetCDF::NetCDF(metFile_idim, metFile_jdim, metFile_kdim) + { + NLON = metFile_idim; + NLAT = metFile_jdim; + NALT = metFile_kdim; + + longitude = new float [NLON]; + latitude = new float [NLAT]; + altitude = new float [NALT]; + u = new float[NALT*NLON*NLAT]; + v = new float[NALT*NLON*NLAT]; + w = new float[NALT*NLON*NLAT]; + dudx = new float[NALT*NLON*NLAT]; + dvdx = new float[NALT*NLON*NLAT]; + dwdx = new float[NALT*NLON*NLAT]; + dudy = new float[NALT*NLON*NLAT]; + dvdy = new float[NALT*NLON*NLAT]; + dwdy = new float[NALT*NLON*NLAT]; + dudz = new float[NALT*NLON*NLAT]; + dvdz = new float[NALT*NLON*NLAT]; + dwdz = new float[NALT*NLON*NLAT]; + thetarhobar = new float[NALT*NLON*NLAT]; + dpibardx = new float[NALT*NLON*NLAT]; + dpibardy = new float[NALT*NLON*NLAT]; + + } + +NetCDF_XYZ::~NetCDF_XYZ() +{ + delete[] longitude; + delete[] latitude; + delete[] altitude; + delete[] u; + delete[] v; + delete[] w; + delete[] dudx; + delete[] dvdx; + delete[] dwdx; + delete[] dudy; + delete[] dvdy; + delete[] dwdy; + delete[] dudz; + delete[] dvdz; + delete[] dwdz; + delete[] thetarhobar; + delete[] dpibardx; + delete[] dpibardy; + +} + +int NetCDF_XYZ::readNetCDF(const char* filename) { + NcError err(NcError::verbose_nonfatal); + NcFile dataFile(filename, NcFile::ReadOnly); + + if(!dataFile.is_valid()) + return NC_ERR; + + // Get pointers to the latitude and longitude variables. + NcVar *lonVar, *latVar, *altVar; + if (!(lonVar = dataFile.get_var("lon"))) + return NC_ERR; + if (!(latVar = dataFile.get_var("lat"))) + return NC_ERR; + if (!(altVar = dataFile.get_var("z"))) + return NC_ERR; + + // Get the lat/lon data from the file. + if (!lonVar->get(longitude, NLON)) + return NC_ERR; + if (!latVar->get(latitude, NLAT)) + return NC_ERR; + if (!altVar->get(altitude, NALT)) + return NC_ERR; + + // Get pointers to the pressure and temperature variables. + NcVar *uVar,*vVar,*wVar,*dudxVar,*dvdxVar,*dwdxVar,*dudyVar,*dvdyVar,*dwdyVar,*dudzVar,*dvdzVar,*dwdzVar, *trbVar, *dpibdxVar, *dpibdyVar; + + if (!(uVar = dataFile.get_var("u"))) + return NC_ERR; + if (!(vVar = dataFile.get_var("v"))) + return NC_ERR; + if (!(wVar = dataFile.get_var("w"))) + return NC_ERR; + if (!(dudxVar = dataFile.get_var("dudx"))) + return NC_ERR; + if (!(dvdxVar = dataFile.get_var("dvdx"))) + return NC_ERR; + if (!(dwdxVar = dataFile.get_var("dwdx"))) + return NC_ERR; + if (!(dudyVar = dataFile.get_var("dudy"))) + return NC_ERR; + if (!(dvdyVar = dataFile.get_var("dvdy"))) + return NC_ERR; + if (!(dwdyVar = dataFile.get_var("dwdy"))) + return NC_ERR; + if (!(dudzVar = dataFile.get_var("dudz"))) + return NC_ERR; + if (!(dvdzVar = dataFile.get_var("dvdz"))) + return NC_ERR; + if (!(dwdzVar = dataFile.get_var("dwdz"))) + return NC_ERR; + if (!(trbVar = dataFile.get_var("trb"))) + return NC_ERR; + if (!(dpibdxVar = dataFile.get_var("dpibdx"))) + return NC_ERR; + if (!(dpibdyVar = dataFile.get_var("dpibdy"))) + return NC_ERR; + + if (!uVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!vVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!wVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dudxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dvdxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dwdxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dudyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dvdyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dwdyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dudzVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dvdzVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dwdzVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!trbVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dpibdxVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + if (!dpibdyVar->set_cur(NREC, 0, 0, 0)) + return NC_ERR; + + if (!uVar->get(u, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!vVar->get(v, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!wVar->get(w, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dudxVar->get(dudx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dvdxVar->get(dvdx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dwdxVar->get(dwdx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dudyVar->get(dudy, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dvdyVar->get(dvdy, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dwdyVar->get(dwdy, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dudzVar->get(dudz, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dvdzVar->get(dvdz, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dwdzVar->get(dwdz, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!trbVar->get(thetarhobar, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dpibdxVar->get(dpibardx, 1, NALT, NLAT, NLON)) + return NC_ERR; + if (!dpibdyVar->get(dpibardy, 1, NALT, NLAT, NLON)) + return NC_ERR; + + return 0; +} + + +double NetCDF_XYZ::getValue(const int &i,const int &j,const int &k, const QString &varName) +{ + //Returned values are all in SI units + double value_out; + + if ((i<0) or (i > NLON-1)) return false; + if ((j<0) or (j > NLAT-1)) return false; + if ((k<0) or (k > NALT-1)) return false; + + if (varName=="lon") { + value_out= longitude[i]; + } else if (varName=="lat") { + value_out= latitude[j]; + } else if (varName=="z") { + value_out= altitude[k]*1000.0; + } else if (varName=="u") { + value_out= u[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="v") { + value_out= v[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="w") { + value_out= w[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="dudx") { + value_out= dudx[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dvdx") { + value_out= dvdx[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dwdx") { + value_out= dwdx[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dudy") { + value_out= dudy[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dvdy") { + value_out= dvdy[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dwdy") { + value_out= dwdy[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dudz") { + value_out= dudz[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dvdz") { + value_out= dvdz[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="dwdz") { + value_out= dwdz[k*NLAT*NLON + j*NLAT + i]/1.0E5; + } else if (varName=="trb") { + value_out= thetarhobar[k*NLAT*NLON + j*NLAT + i]; + } else if (varName=="dpibdx") { + value_out= dpibardx[k*NLAT*NLON + j*NLAT + i]/1000.0; + } else if (varName=="dpibdy") { + value_out= dpibardy[k*NLAT*NLON + j*NLAT + i]/1000.0; + } else if (varName=="A") { + value_out= this->calc_A(i,j,k); + } else if (varName=="B") { + value_out= this->calc_B(i,j,k); + } else if (varName=="C") { + value_out= this->calc_C(i,j,k); + + } else { + std::cout << "Requested Variable unknown. \n"; + std::cout << varName.toStdString() << "\n"; + return false; + } + + return value_out; +} + + +double NetCDF_XYZ::calc_A(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,(QString)"trb"); + double u = this->getValue(i,j,k,(QString)"u"); + double dudx = this->getValue(i,j,k,(QString)"dudx"); + double v = this->getValue(i,j,k,(QString)"v"); + double dudy = this->getValue(i,j,k,(QString)"dudy"); + double w = this->getValue(i,j,k,(QString)"w"); + double dudz = this->getValue(i,j,k,(QString)"dudz"); + double dpibdx = this->getValue(i,j,k,(QString)"dpibdx"); + float c_p = 1005.7; + + if (thetarhobar==-999 or u==-999 or dudx*1.0E5==-999 or v==-999 or dudy*1.0E5==-999 or w==-999 or dudz*1.0E5==-999 or dpibdx*1000.0==-999){ + return -999;} + + double a = 1.0/(c_p*thetarhobar)* (u*dudx+v*dudy+w*dudz-f*v)+dpibdx; + return a; +} + + +double NetCDF_XYZ::calc_B(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,(QString)"trb"); + double u = this->getValue(i,j,k,(QString)"u"); + double dvdx = this->getValue(i,j,k,(QString)"dvdx"); + double v = this->getValue(i,j,k,(QString)"v"); + double dvdy = this->getValue(i,j,k,(QString)"dvdy"); + double w = this->getValue(i,j,k,(QString)"w"); + double dvdz = this->getValue(i,j,k,(QString)"dvdz"); + double dpibdy = this->getValue(i,j,k,(QString)"dpibdy"); + float c_p = 1005.7; + + if (thetarhobar==-999 or u==-999 or dvdx*1.0E5==-999 or v==-999 or dvdy*1.0E5==-999 or w==-999 or dvdz*1.0E5==-999 or dpibdy*1000.0==-999){ + return -999;} + + double b = 1.0/(c_p*thetarhobar)*(u*dvdx+v*dvdy+w*dvdz+f*u)+dpibdy; // trp neglected + return b; +} + +double NetCDF_XYZ::calc_C(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,(QString)"trb"); + double u = this->getValue(i,j,k,(QString)"u"); + double dwdx = this->getValue(i,j,k,(QString)"dwdx"); + double v = this->getValue(i,j,k,(QString)"v"); + double dwdy = this->getValue(i,j,k,(QString)"dwdy"); + double w = this->getValue(i,j,k,(QString)"w"); + double dwdz = this->getValue(i,j,k,(QString)"dwdz"); + float c_p = 1005.7; + float g = 9.81; + + if (thetarhobar==-999 or u==-999 or dwdx*1.0E5==-999 or v==-999 or dwdy*1.0E5==-999 or w==-999 or dwdz*1.0E5==-999){ + return -999;} + + double c = 1.0/(c_p*thetarhobar)*(u*dwdx+v*dwdy+w*dwdz); + return c; +} + +double NetCDF_XYZ::calc_D(const int &i,const int &j,const int &k) +{ + + double thetarhobar = this->getValue(i,j,k,(QString)"trb"); + double dAdz = this->getDerivative(i,j,k,(QString)"A",3)*1000.0; // per km + double dCdx = this->getDerivative(i,j,k,(QString)"C",1)*1000.0; // per km + float c_p = 1005.7; + float g = 9.81; + + if (thetarhobar==-999 or dAdz==-999000 or dCdx==-999000){ + return -999;} + + double d = (dAdz-dCdx)*thetarhobar*thetarhobar*(-c_p/g); + + return d; +} + +double NetCDF_XYZ::calc_E(const int &i,const int &j,const int &k) +{ + double thetarhobar = this->getValue(i,j,k,(QString)"trb"); + double dBdz = this->getDerivative(i,j,k,(QString)"B",3)*1000.0; // per km + double dCdy = this->getDerivative(i,j,k,(QString)"C",2)*1000.0; // per km + float c_p = 1005.7; + float g = 9.81; + + if (thetarhobar==-999 or dBdz==-999000 or dCdy==-999000){ + return -999;} + + double e = (dBdz-dCdy)*thetarhobar*thetarhobar*(-c_p/g); + + return e; +} + + +double NetCDF_XYZ::getDerivative(const int &i,const int &j,const int &k, const QString &var, const int &der) +{ + // input variable "der" specifies direction of derivation: 1=dx, 2=dy, 3=dz + double derivative; + QString derDir; + // Geographic functions + GeographicLib::TransverseMercatorExact tm = GeographicLib::TransverseMercatorExact::UTM(); + double referenceLon = -90.0; //arbitrary + double x1,x2,y1,y2; + + switch ( der ) { + case 1: + derDir = "X"; + if (i==0){ + //need to convert lat/lon differences to kms + tm.Forward(referenceLon,this->getValue(i+1,j,k,"lat"),this->getValue(i+1,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i+1,j,k,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i+1,j,k,var)-this->getValue(i,j,k,var))/distance; + } + } else if (i== NLON-1) { + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i-1,j,k,"lat"),this->getValue(i-1,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j,k,var)==-999 || this->getValue(i-1,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k,var)-this->getValue(i-1,j,k,var))/distance; + } + } else { + tm.Forward(referenceLon,this->getValue(i+1,j,k,"lat"),this->getValue(i+1,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i-1,j,k,"lat"),this->getValue(i-1,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i+1,j,k,var)==-999 || this->getValue(i-1,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i+1,j,k,var)-this->getValue(i-1,j,k,var))/distance; + } + } + break; + case 2: + derDir = "Y"; + if (j==0){ + tm.Forward(referenceLon,this->getValue(i,j+1,k,"lat"),this->getValue(i,j+1,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j+1,k,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j+1,k,var)-this->getValue(i,j,k,var))/distance; + } + } else if (j== NLAT-1) { + tm.Forward(referenceLon,this->getValue(i,j,k,"lat"),this->getValue(i,j,k,"lon"),x1,y1); + tm.Forward(referenceLon,this->getValue(i,j-1,k,"lat"),this->getValue(i,j-1,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j,k,var)==-999 || this->getValue(i,j-1,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k,var)-this->getValue(i,j-1,k,var))/distance; + } + } else { + tm.Forward(referenceLon,this->getValue(i,j+1,k,"lat"),this->getValue(i,j+1,k,"lon"),x2,y2); + tm.Forward(referenceLon,this->getValue(i,j-1,k,"lat"),this->getValue(i,j-1,k,"lon"),x2,y2); + double distance = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); + if (this->getValue(i,j+1,k,var)==-999 || this->getValue(i,j-1,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j+1,k,var)-this->getValue(i,j-1,k,var))/distance; + } + } + break; + case 3: + derDir = "Z"; + if (k==0){ + if (this->getValue(i,j,k+1,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k+1,var)-this->getValue(i,j,k,var))/(this->getValue(i,j,k+1,"z")-this->getValue(i,j,k,"z")); + } + } else if (k== NALT-1) { + if (this->getValue(i+1,j,k,var)==-999 || this->getValue(i,j,k,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k,var)-this->getValue(i,j,k-1,var))/(this->getValue(i,j,k,"z")-this->getValue(i,j,k-1,"z")); + } + } else { + if (this->getValue(i,j,k+1,var)==-999 || this->getValue(i,j,k-1,var)==-999){ + derivative = -999; + } else { + derivative = (this->getValue(i,j,k+1,var)-this->getValue(i,j,k-1,var))/(this->getValue(i,j,k+1,"z")-this->getValue(i,j,k-1,"z")); + } + } + break; + default: + std::cout << "Unknown value for calculating derivative. Valid options are 1, 2 and 3.\n"; + exit(1); + } + return derivative; +} diff --git a/src/thermo_reference/NetCDF_XYZ.h b/src/thermo_reference/NetCDF_XYZ.h new file mode 100644 index 0000000..749b193 --- /dev/null +++ b/src/thermo_reference/NetCDF_XYZ.h @@ -0,0 +1,62 @@ + /* + * NetCDF.h + * samurai + * + * Created by Annette Foerster. + * Based on example code from netCDF library + * Copyright 2013 Michael Bell. All rights reserved. + * + */ + +#ifndef NETCDF_XYZ_H +#define NETCDF_XYZ_H + +#include "NetCDF.h" +#include +#include + +class NetCDF_XYZ : public NetCDF +{ + +public: + NetCDF_XYZ(const int& metFile_idim, const int& metFile_jdim, const int& metFile_kdim); + + ~NetCDF_XYZ(); + + int readNetCDF(const char* filename); + double getValue(const int &i,const int &j,const int &k,const QString &varName); + double getDerivative(const int &i,const int &j,const int &k, const QString &var, const int &der); + double calc_A(const int &i,const int &j,const int &k); + double calc_B(const int &i,const int &j,const int &k); + double calc_C(const int &i,const int &j,const int &k); + double calc_D(const int &i,const int &j,const int &k); + double calc_E(const int &i,const int &j,const int &k); + +protected: + int NLON, NLAT; + + float* longitude; + float* latitude; + float* altitude; + float* u; + float* v; + float* w; + float* dudx; + float* dvdx; + float* dwdx; + float* dudy; + float* dvdy; + float* dwdy; + float* dudz; + float* dvdz; + float* dwdz; + float* thetarhobar; + float* dpibardx; + float* dpibardy; + + QString varName; + +}; + + +#endif diff --git a/src/thermo_reference/VarDriverThermo.cpp b/src/thermo_reference/VarDriverThermo.cpp new file mode 100644 index 0000000..fe24fc3 --- /dev/null +++ b/src/thermo_reference/VarDriverThermo.cpp @@ -0,0 +1,624 @@ +/* + * VarDriverThermo.cpp + * samurai + * + * Copyright 2008 Michael Bell. All rights reserved. + * + */ + +#include "VarDriverThermo.h" +#include +#include +#include +#include +#include +#include +#include +#include "RecursiveFilter.h" +#include + +// Constructor +VarDriverThermo::VarDriverThermo() +: VarDriver3D() +{ + numVars = 7; + numDerivatives = 4; + obMetaSize = 7; +} + +// Destructor +VarDriverThermo::~VarDriverThermo() +{ +} + +/* This routine is the main initializer of the analysis */ + +bool VarDriverThermo::initialize(const QDomElement& configuration) +{ + // Run a 3D vortex background field + cout << "Initializing SAMURAI Thermo" << endl; + + // Parse the XML configuration file + if (!parseXMLconfig(configuration)) return false; + + // Validate the 3D specific parameters + if (!validateXMLconfig()) return false; + + // Validate the run geometry + if (configHash.value("mode") == "XYZ") { + runMode = XYZ; + } else if (configHash.value("mode") == "RTZ") { + runMode = RTZ; + } else { + cout << "Unrecognized run mode " << configHash.value("mode").toStdString() << ", Aborting...\n"; + return false; + } + + // Define the grid dimensions + imin = configHash.value("i_min").toFloat(); + imax = configHash.value("i_max").toFloat(); + iincr = configHash.value("i_incr").toFloat(); + idim = (int)((imax - imin)/iincr) + 1; + + jmin = configHash.value("j_min").toFloat(); + jmax = configHash.value("j_max").toFloat(); + jincr = configHash.value("j_incr").toFloat(); + jdim = (int)((jmax - jmin)/jincr) + 1; + + kmin = configHash.value("k_min").toFloat(); + kmax = configHash.value("k_max").toFloat(); + kincr = configHash.value("k_incr").toFloat(); + kdim = (int)((kmax - kmin)/kincr) + 1; + + // The recursive filter uses a fourth order stencil to spread the observations, 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"; + return false; + } + if (jdim < 4) { + 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"; + 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"; + + int uStateSize = 8*(idim+1)*(jdim+1)*(kdim+1)*(numVars); + int bStateSize = (idim+2)*(jdim+2)*(kdim+2)*numVars; + cout << "Physical (mish) State size = " << uStateSize << "\n"; + cout << "Nodal State size = " << bStateSize << ", Grid dimensions:\n"; + + // Load the BG into a empty vector + bgU = new real[uStateSize]; + bgWeights = new real[uStateSize]; + for (int i=0; i < uStateSize; i++) { + bgU[i] = 0.0; + bgWeights[i] = 0.0; + } + + /* Set the data path */ + dataPath.setPath(configHash.value("data_directory")); + if (dataPath.exists()) { + dataPath.setFilter(QDir::Files); + dataPath.setSorting(QDir::Name); + } else { + cout << "Can't find data directory: " << configHash.value("data_directory").toStdString() << endl; + return false; + } + + /* Check to make sure the output path exists */ + QDir outputPath(configHash.value("output_directory")); + if (!outputPath.exists()) { + cout << "Can't find output directory: " << configHash.value("output_directory").toStdString() << endl; + return false; + } + + // Define the Reference state + QString refSounding = dataPath.absoluteFilePath(configHash.value("ref_state")); + refstate = new ReferenceState(refSounding); + cout << "Reference profile: Z\t\tQv\tRhoa\tRho\tH\tTemp\tPressure\n"; + for (real k = kmin; k < kmax+kincr; k+= kincr) { + cout << " " << k << "\t"; + for (int i = 0; i < 6; i++) { + //real var = refstate->getReferenceVariable(i, k*1000); + //if (i == 0) var = refstate->bhypInvTransform(var); + real var = 0.0; + cout << setw(9) << setprecision(4) << var << "\t"; + } + cout << "\n"; + } + cout << setprecision(9); + + // Read in the Frame centers + // Ideally, create a time-based spline from limited center fixes here + // but just load 1 second centers into vector for now + readFrameCenters(); + + // Get the reference center + QTime reftime = QTime::fromString(configHash.value("ref_time"), "hh:mm:ss"); + QString refstring = reftime.toString(); + bool foundref = false; + for (unsigned int fi = 0; fi < frameVector.size(); fi++) { + QDateTime frametime = frameVector[fi].getTime(); + if (reftime == frametime.time()) { + QString tempstring; + QDate refdate = frametime.date(); + QDateTime unixtime(refdate, reftime, Qt::UTC); + configHash.insert("ref_lat", tempstring.setNum(frameVector[fi].getLat())); + configHash.insert("ref_lon", tempstring.setNum(frameVector[fi].getLon())); + configHash.insert("ref_time", tempstring.setNum(unixtime.toTime_t())); + cout << "Found matching reference time " << refstring.toStdString() + << " at " << frameVector[fi].getLat() << ", " << frameVector[fi].getLon() << "\n"; + foundref = true; + break; + } + } + if (!foundref) { + cout << "Error finding reference time, please check date and time in XML file\n"; + return false; + } + + /* Set the maximum number of iterations to the multipass reduction factor + Multiple outer loops will reduce the cutoff wavelengths and background error variance */ + maxIter = configHash.value("num_iterations").toInt(); + + /* Optionally load a set of background estimates and interpolate to the Gaussian mish */ + + metFile = configHash.value("infile"); + metFile_idim = configHash.value("infile_idim").toInt(); + metFile_jdim = configHash.value("infile_jdim").toInt(); + metFile_kdim = configHash.value("infile_kdim").toInt(); + + if(!this->loadObservations(metFile,metFile_idim,metFile_jdim,metFile_kdim, &obVector)) { + // For testing purposes, comment out line above and use this one instead: if(!this->testing(&obVector)) { + cout << "Loading Observations failed ...Exit." << endl; + return EXIT_FAILURE; + } + + cout << "Number of New Observations: " << obVector.size() << endl; + + + // 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); + + } + } + + } + + +// AF for now set the background to zero and don't allow any additional observations, i.e. skip loadBackgroundObs, adjustBackground, preProcessMetObs and loadMetObs + + // We are done with the bgWeights, so free up that memory + delete[] bgWeights; + + if (runMode == XYZ) { + obCost3D = new CostFunctionThermoXYZ(obVector.size(), bStateSize); + } else if (runMode == RTZ) { + obCost3D = new CostFunctionThermoRTZ(obVector.size(), bStateSize); + } + obCost3D->initialize(&configHash, bgU, obs, refstate); + + + // If we got here, then everything probably went OK! + return true; +} + +/* This routine drives the CostFunction minimization + There is support for an outer loop to change the background + error covariance or update non-linear observation operators */ + +bool VarDriverThermo::run() +{ + int iter=1; + while (iter <= maxIter) { + cout << "Outer Loop Iteration: " << iter << endl; + obCost3D->initState(iter); + obCost3D->minimize(); + obCost3D->updateBG(); + iter++; + + // Optionally update the analysis parameters for an additional iteration + updateAnalysisParams(iter); + } + + return true; + +} + +/* Clean up all that allocated memory */ + +bool VarDriverThermo::finalize() +{ + obCost3D->finalize(); + delete[] obs; + delete[] bgU; + delete obCost3D; + delete refstate; + return true; +} + +bool VarDriverThermo::loadObservations(QString& metFile,const int &metFile_idim,const int &metFile_jdim,const int &metFile_kdim , QList* obVector) +{ + 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); + } + + cout << "Read in NetCDF File ... " << endl; + if (ncFile->readNetCDF(metFile.toAscii().data()) != 0) { + cout << "Error reading NetCDF file\n"; + exit(1); + } + + cout << "Load Observations ... " << endl; + + int nalt = metFile_kdim; + int nx = metFile_idim; + int ny = metFile_jdim; + + 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); + if ((fi < 0) or (fi > (int)frameVector.size())) { + cout << "Time problem with observation " << fi << endl; + exit(1); + } + + if (configHash.value("mode") == "RTZ") { + cout << "RTZ run mode not implemented yet " << endl; + + } else if (configHash.value("mode") == "XYZ") { + std::cout << "Entered XYZ run mode " << endl; + + int nlon = nx; + int nlat = ny; + 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,(QString)"lon"); + double lat = ncFile->getValue(i,j,k,(QString)"lat"); + double alt = ncFile->getValue(i,j,k,(QString)"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 = configHash.value("ref_lon").toFloat(); + 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(101); + varOb.setCartesianX(obX); + varOb.setCartesianY(obY); + varOb.setAltitude(alt_km); + varOb.setTime(obTime.toTime_t()); + + + // 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,(QString)"trb"); + + + double u = ncFile->getValue(i,j,k,(QString)"u"); + double v = ncFile->getValue(i,j,k,(QString)"v"); + double w = ncFile->getValue(i,j,k,(QString)"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){ + 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(configHash.value("thermo_A_error").toFloat()); + obVector->push_back(varOb); + varOb.setWeight(0,0,1); + + varOb.setOb(b*1E3*1E3); + varOb.setWeight(-0.001*1E3,0,2); + varOb.setError(configHash.value("thermo_B_error").toFloat()); + 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(configHash.value("thermo_C_error").toFloat()); + obVector->push_back(varOb); + varOb.setWeight(0,0,3); + varOb.setWeight(0,1,0); + + varOb.setOb(d); + varOb.setWeight(1,1,1); + varOb.setError(configHash.value("thermo_D_error").toFloat()); + obVector->push_back(varOb); + varOb.setWeight(0,1,1); + + varOb.setOb(e); + varOb.setWeight(1,1,2); + varOb.setError(configHash.value("thermo_E_error").toFloat()); + obVector->push_back(varOb); + varOb.setWeight(0,1,2); + + } + } + } + + } else { + cout << "Unrecognized run mode " << configHash.value("mode").toStdString() << ", Aborting...\n"; + return false; + } + + return true; + +} + +bool VarDriverThermo::validateXMLconfig() +{ + + // Validate the hash -- multiple passes are not validated currently + QStringList configKeys; + configKeys << "i_min" << "i_max" << "i_incr" << + "j_min" << "j_max" << "j_incr" << + "k_min" << "k_max" << "k_incr" << + "i_pip_bcL" << "i_pip_bcR" << "j_pip_bcL" << "j_pip_bcR" << "k_pip_bcL" << "k_pip_bcR" << + "i_thetarhop_bcL" << "i_thetarhop_bcR" << "j_thetarhop_bcL" << "j_thetarhop_bcR" << "k_thetarhop_bcL" << "k_thetarhop_bcR" << + "data_directory" << "output_directory"; + for (int i = 0; i < configKeys.count(); i++) { + if (!configHash.contains(configKeys.at(i))) { + cout << "No configuration found for <" << configKeys.at(i).toStdString() << "> aborting..." << endl; + return false; + } + } + return true; +} + +bool VarDriverThermo::testing(QList* obVector) +{ + // 3 x_1 + 2 x_2 - 1 x_3 = 1 + // 2 x_1 - 2 x_2 + 4 x_3 = -2 + //-1 x_1 + 0.5 x_2 & -1 x_3 = 0 + + //Solution: x_1 = 1, x_2 = -2, x_3 = -2 + + QString file,datestr,timestr; + file = "20050921_18.nc"; + 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 { + 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); + if ((fi < 0) or (fi > (int)frameVector.size())) { + cout << "Time problem with observation " << fi << endl; + exit(1); + } + +//for (int i = -10; i < 10; ++i) { + //for (int j = -10; j < 10; ++j) { + //for (int k = 0; k < 7; ++k) { + + int i=0; + int j=0; + int k=3; + + Observation varOb; + + varOb.setType(101); + varOb.setCartesianX(i); + varOb.setCartesianY(j); + varOb.setAltitude(k); + varOb.setTime(obTime.toTime_t()); + + + // 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); + } + } + + varOb.setOb(1); + varOb.setWeight(3,0,0); + varOb.setWeight(2,1,0); + varOb.setWeight(-1,2,0); + varOb.setError(0.1); + obVector->push_back(varOb); + varOb.setWeight(0,0,0); + varOb.setWeight(0,1,0); + varOb.setWeight(0,2,0); + + + varOb.setOb(-2); + varOb.setWeight(2,0,0); + varOb.setWeight(-2,1,0); + varOb.setWeight(4,2,0); + varOb.setError(0.1); + obVector->push_back(varOb); + varOb.setWeight(0,0,0); + varOb.setWeight(0,1,0); + varOb.setWeight(0,2,0); + + varOb.setOb(0); + varOb.setWeight(-1,0,0); + varOb.setWeight(0.5,1,0); + varOb.setWeight(-1,2,0); + varOb.setError(0.1); + obVector->push_back(varOb); + varOb.setWeight(0,0,0); + varOb.setWeight(0,1,0); + varOb.setWeight(0,2,0); + + + // } + //} +//} + + return true; +} + +bool VarDriverThermo::testing_rtz(QList* obVector) +{ + + QString file,datestr,timestr; + file = "20050921_18.nc"; + 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 { + 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); + if ((fi < 0) or (fi > (int)frameVector.size())) { + cout << "Time problem with observation " << fi << endl; + exit(1); + } + + int i=5; + int j=90; + int k=3; + + Observation varOb; + + varOb.setType(101); + varOb.setRadius(i); + varOb.setTheta(j); + varOb.setAltitude(k); + varOb.setTime(obTime.toTime_t()); + + // 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); + } + } + + varOb.setOb(1); + varOb.setWeight(1,1,0); + varOb.setError(0.0001); + obVector->push_back(varOb); + varOb.setWeight(0,0,0); + + + /*varOb.setOb(2); + varOb.setWeight(1,1,0); + varOb.setError(0.01); + obVector->push_back(varOb); + varOb.setWeight(0,1,0); + + varOb.setOb(3); + varOb.setWeight(1,2,0); + varOb.setError(0.01); + obVector->push_back(varOb); + varOb.setWeight(0,2,0);*/ + + return true; +} diff --git a/src/thermo_reference/VarDriverThermo.h b/src/thermo_reference/VarDriverThermo.h new file mode 100644 index 0000000..a7991ad --- /dev/null +++ b/src/thermo_reference/VarDriverThermo.h @@ -0,0 +1,65 @@ +/* + * VarDriverThermo.h + * samurai + * + * Created by Annette Foerster on 8/7/13. + * Copyright 2008 Michael Bell. All rights reserved. + * + */ + +#ifndef VARDRIVERTHERMO_H +#define VARDRIVERTHERMO_H + +#include "VarDriver3D.h" +#include "VarDriver.h" +#include "BSpline.h" +#include "Observation.h" +#include "CostFunctionThermoXYZ.h" +#include "CostFunctionThermoRTZ.h" +#include "CostFunctionThermo.h" +#include "MetObs.h" +#include "FrameCenter.h" +#include "NetCDF.h" +#include "NetCDF_RTZ.h" +#include "NetCDF_XYZ.h" +#include +#include +#include +#include +#include +#include + +using namespace std; + +class VarDriverThermo : public VarDriver3D +{ + +public: + VarDriverThermo(); + ~VarDriverThermo(); + + // ESMF type calls + bool initialize(const QDomElement& configuration); + bool run(); + bool finalize(); + bool validateXMLconfig(); + bool testing(QList* obVector); + bool testing_rtz(QList* obVector); + + + +private: + int test; + CostFunctionThermo* obCost3D; + bool loadObservations(QString& metFile,const int &metFile_idim,const int &metFile_jdim,const int &metFile_kdim , QList* obVector); + NetCDF* ncFile; + QString metFile; + int metFile_idim; + int metFile_jdim; + int metFile_kdim; + QList obVector; + + real* obs; +}; + +#endif