test_htgop_load.c (17005B)
1 /* Copyright (C) 2018-2021, 2023 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #define _POSIX_C_SOURCE 200112L /* nextafter support */ 17 18 #include "htgop.h" 19 #include "htgop_reader.h" 20 21 #include "test_htgop_utils.h" 22 23 #include <rsys/cstr.h> 24 25 static void 26 check_position_to_layer_id(const struct htgop* htgop) 27 { 28 struct htgop_level lvl0, lvl1; 29 const size_t N = 10000; 30 double org, sz; 31 size_t ilay; 32 size_t i, n; 33 34 CHK(htgop_get_levels_count(htgop, &n) == RES_OK); 35 CHK(n > 1); 36 37 /* Test the htgop_position_to_layer_id function */ 38 CHK(htgop_get_level(htgop, 0, &lvl0) == RES_OK); 39 CHK(htgop_get_level(htgop, n-1, &lvl1) == RES_OK); 40 org = lvl0.height; 41 sz = lvl1.height - lvl0.height; 42 43 CHK(htgop_position_to_layer_id(NULL, lvl0.height, &ilay) == RES_BAD_ARG); 44 CHK(htgop_position_to_layer_id 45 (htgop, nextafter(lvl0.height,-DBL_MAX), &ilay) == RES_BAD_ARG); 46 CHK(htgop_position_to_layer_id 47 (htgop, nextafter(lvl1.height, DBL_MAX), &ilay) == RES_BAD_ARG); 48 CHK(htgop_position_to_layer_id(htgop, lvl0.height, NULL) == RES_BAD_ARG); 49 50 CHK(htgop_position_to_layer_id(htgop, lvl0.height, &ilay) == RES_OK); 51 CHK(ilay == 0); 52 CHK(htgop_position_to_layer_id(htgop, lvl1.height, &ilay) == RES_OK); 53 CHK(ilay == n - 2); 54 55 FOR_EACH(i, 0, N) { 56 const double r = rand_canonic(); 57 const double pos = org + r * sz; 58 59 CHK(htgop_position_to_layer_id(htgop, pos, &ilay) == RES_OK); 60 CHK(htgop_get_level(htgop, ilay+0, &lvl0) == RES_OK); 61 CHK(htgop_get_level(htgop, ilay+1, &lvl1) == RES_OK); 62 CHK(pos >= lvl0.height); 63 CHK(pos <= lvl1.height); 64 } 65 } 66 67 #define DOMAIN sw 68 #include "test_htgop_check_specints.h" 69 70 #define DOMAIN lw 71 #include "test_htgop_check_specints.h" 72 73 int 74 main(int argc, char** argv) 75 { 76 struct mem_allocator allocator; 77 struct reader rdr; 78 struct htgop* htgop; 79 struct htgop_ground ground; 80 struct htgop_layer lay; 81 struct htgop_level lvl; 82 struct htgop_spectral_interval specint; 83 struct htgop_layer_lw_spectral_interval lw_specint; 84 struct htgop_layer_lw_spectral_interval_tab lw_tab; 85 struct htgop_layer_sw_spectral_interval sw_specint; 86 struct htgop_layer_sw_spectral_interval_tab sw_tab; 87 size_t nlays, nlvls, nspecints_lw, nspecints_sw; 88 size_t ilay, ilvl, itab, iquad, ispecint; 89 size_t tab_len; 90 const double* wnums; 91 FILE* fp; 92 unsigned long ul; 93 double dbl; 94 95 if(argc < 2) { 96 fprintf(stderr, "Usage: %s FILENAME\n", argv[0]); 97 return 1; 98 } 99 100 CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); 101 102 CHK((fp = fopen(argv[1], "r")) != NULL); 103 reader_init(&rdr, fp, argv[1]); 104 105 CHK(htgop_create(NULL, &allocator, 1, &htgop) == RES_OK); 106 CHK(htgop_load(NULL, NULL) == RES_BAD_ARG); 107 CHK(htgop_load(htgop, NULL) == RES_BAD_ARG); 108 CHK(htgop_load(NULL, argv[1]) == RES_BAD_ARG); 109 CHK(htgop_load(htgop, argv[1]) == RES_OK); 110 111 /* #levels */ 112 CHK(htgop_get_levels_count(NULL, NULL) == RES_BAD_ARG); 113 CHK(htgop_get_levels_count(htgop, NULL) == RES_BAD_ARG); 114 CHK(htgop_get_levels_count(NULL, &nlvls) == RES_BAD_ARG); 115 CHK(htgop_get_levels_count(htgop, &nlvls) == RES_OK); 116 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 117 CHK(nlvls == ul); 118 119 /* #layers */ 120 CHK(htgop_get_layers_count(NULL, NULL) == RES_BAD_ARG); 121 CHK(htgop_get_layers_count(htgop, NULL) == RES_BAD_ARG); 122 CHK(htgop_get_layers_count(NULL, &nlays) == RES_BAD_ARG); 123 CHK(htgop_get_layers_count(htgop, &nlays) == RES_OK); 124 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 125 CHK(nlays == ul); 126 127 /* Ground temperature */ 128 CHK(htgop_get_ground(NULL, NULL) == RES_BAD_ARG); 129 CHK(htgop_get_ground(htgop, NULL) == RES_BAD_ARG); 130 CHK(htgop_get_ground(NULL, &ground) == RES_BAD_ARG); 131 CHK(htgop_get_ground(htgop, &ground) == RES_OK); 132 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 133 CHK(ground.temperature == dbl); 134 135 /* Per level pressure */ 136 CHK(htgop_get_level(NULL, nlvls, NULL) == RES_BAD_ARG); 137 CHK(htgop_get_level(htgop, nlvls, NULL) == RES_BAD_ARG); 138 CHK(htgop_get_level(NULL, 0, NULL) == RES_BAD_ARG); 139 CHK(htgop_get_level(htgop, 0, NULL) == RES_BAD_ARG); 140 CHK(htgop_get_level(NULL, nlvls, &lvl) == RES_BAD_ARG); 141 CHK(htgop_get_level(htgop, nlvls, &lvl) == RES_BAD_ARG); 142 CHK(htgop_get_level(NULL, 0, &lvl) == RES_BAD_ARG); 143 FOR_EACH(ilvl, 0, nlvls) { 144 CHK(htgop_get_level(htgop, ilvl, &lvl) == RES_OK); 145 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 146 CHK(dbl == lvl.pressure); 147 } 148 /* Per level temperature */ 149 FOR_EACH(ilvl, 0, nlvls) { 150 CHK(htgop_get_level(htgop, ilvl, &lvl) == RES_OK); 151 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 152 CHK(dbl == lvl.temperature); 153 } 154 /* Per level height */ 155 FOR_EACH(ilvl, 0, nlvls) { 156 CHK(htgop_get_level(htgop, ilvl, &lvl) == RES_OK); 157 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 158 CHK(dbl == lvl.height); 159 if(ilvl) { /* Check strict ascending order */ 160 CHK(htgop_get_level(htgop, ilvl-1, &lvl) == RES_OK); 161 CHK(lvl.height < dbl); 162 } 163 } 164 165 /* Per layer nominal xH2O */ 166 CHK(htgop_get_layer(NULL, nlays, NULL) == RES_BAD_ARG); 167 CHK(htgop_get_layer(htgop, nlays, NULL) == RES_BAD_ARG); 168 CHK(htgop_get_layer(NULL, 0, NULL) == RES_BAD_ARG); 169 CHK(htgop_get_layer(htgop, 0, NULL) == RES_BAD_ARG); 170 CHK(htgop_get_layer(NULL, nlays, &lay) == RES_BAD_ARG); 171 CHK(htgop_get_layer(htgop, nlays, &lay) == RES_BAD_ARG); 172 CHK(htgop_get_layer(NULL, 0, &lay) == RES_BAD_ARG); 173 FOR_EACH(ilay, 0, nlays) { 174 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 175 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 176 CHK(dbl == lay.x_h2o_nominal); 177 } 178 179 /* # tabulated xH2O for each layer */ 180 tab_len = lay.tab_length; 181 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 182 CHK(tab_len == ul); 183 FOR_EACH(itab, 0, tab_len) { 184 FOR_EACH(ilay, 0, nlays) { 185 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 186 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 187 CHK(dbl == lay.x_h2o_tab[itab]); 188 /* Check strict ascending order */ 189 CHK(!itab || lay.x_h2o_tab[itab] > lay.x_h2o_tab[itab-1]); 190 } 191 } 192 193 /* Ground emissivity */ 194 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 195 CHK(dbl == ground.lw_emissivity); 196 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 197 CHK(dbl == ground.sw_emissivity); 198 199 /* #long wave spectral intervals */ 200 CHK(htgop_get_lw_spectral_intervals_count(NULL, NULL) == RES_BAD_ARG); 201 CHK(htgop_get_lw_spectral_intervals_count(htgop, NULL) == RES_BAD_ARG); 202 CHK(htgop_get_lw_spectral_intervals_count(NULL, &nspecints_lw) == RES_BAD_ARG); 203 CHK(htgop_get_lw_spectral_intervals_count(htgop, &nspecints_lw) == RES_OK); 204 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 205 CHK(ul == nspecints_lw); 206 207 /* #short wave spectral intervals */ 208 CHK(htgop_get_sw_spectral_intervals_count(NULL, NULL) == RES_BAD_ARG); 209 CHK(htgop_get_sw_spectral_intervals_count(htgop, NULL) == RES_BAD_ARG); 210 CHK(htgop_get_sw_spectral_intervals_count(NULL, &nspecints_sw) == RES_BAD_ARG); 211 CHK(htgop_get_sw_spectral_intervals_count(htgop, &nspecints_sw) == RES_OK); 212 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 213 CHK(ul == nspecints_sw); 214 215 /* Search for a LW spectral interval */ 216 CHK(htgop_find_lw_spectral_interval_id(NULL, 0, &ispecint) == RES_BAD_ARG); 217 CHK(htgop_find_lw_spectral_interval_id(htgop, 0, NULL) == RES_BAD_ARG); 218 CHK(htgop_find_lw_spectral_interval_id(htgop, 0, &ispecint) == RES_OK); 219 CHK(ispecint == SIZE_MAX); 220 221 /* Search for a SW spectral interval */ 222 CHK(htgop_find_sw_spectral_interval_id(NULL, 0, &ispecint) == RES_BAD_ARG); 223 CHK(htgop_find_sw_spectral_interval_id(htgop, 0, NULL) == RES_BAD_ARG); 224 CHK(htgop_find_sw_spectral_interval_id(htgop, 0, &ispecint) == RES_OK); 225 CHK(ispecint == SIZE_MAX); 226 227 /* Per LW spectral interval data */ 228 CHK(htgop_get_lw_spectral_intervals_wave_numbers(NULL, NULL) == RES_BAD_ARG); 229 CHK(htgop_get_lw_spectral_intervals_wave_numbers(htgop, NULL) == RES_BAD_ARG); 230 CHK(htgop_get_lw_spectral_intervals_wave_numbers(NULL, &wnums) == RES_BAD_ARG); 231 CHK(htgop_get_lw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK); 232 CHK(htgop_get_lw_spectral_interval(NULL, nspecints_lw, NULL) == RES_BAD_ARG); 233 CHK(htgop_get_lw_spectral_interval(htgop, nspecints_lw, NULL) == RES_BAD_ARG); 234 CHK(htgop_get_lw_spectral_interval(NULL, 0, NULL) == RES_BAD_ARG); 235 CHK(htgop_get_lw_spectral_interval(htgop, 0, NULL) == RES_BAD_ARG); 236 CHK(htgop_get_lw_spectral_interval(NULL, nspecints_lw, &specint) == RES_BAD_ARG); 237 CHK(htgop_get_lw_spectral_interval(htgop, nspecints_lw, &specint) == RES_BAD_ARG); 238 CHK(htgop_get_lw_spectral_interval(NULL, 0, &specint) == RES_BAD_ARG); 239 FOR_EACH(ispecint, 0, nspecints_lw) { 240 size_t i; 241 CHK(htgop_get_lw_spectral_interval(htgop, ispecint, &specint) == RES_OK); 242 CHK(specint.wave_numbers[0] == wnums[ispecint+0]); 243 CHK(specint.wave_numbers[1] == wnums[ispecint+1]); 244 245 /* Check the finding of a LW spectral interval */ 246 FOR_EACH(i, 0, 10) { 247 double wnum; 248 size_t id; 249 wnum = wnums[ispecint+0] 250 + rand_canonic() * (wnums[ispecint+1] - wnums[ispecint+0]); 251 CHK(htgop_find_lw_spectral_interval_id(htgop, wnum, &id) == RES_OK); 252 CHK(id == ispecint); 253 } 254 255 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 256 CHK(specint.wave_numbers[0] == dbl); 257 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 258 CHK(specint.wave_numbers[1] == dbl); 259 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 260 CHK(specint.quadrature_length == ul); 261 262 FOR_EACH(iquad, 0, specint.quadrature_length) { 263 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 264 CHK(specint.quadrature_pdf[iquad] == dbl); 265 if(iquad == 0) { 266 CHK(specint.quadrature_cdf[iquad] == dbl); 267 } else { 268 const double pdf = 269 specint.quadrature_cdf[iquad] - specint.quadrature_cdf[iquad-1]; 270 CHK(eq_eps(pdf, dbl, 1.e-6)); 271 } 272 } 273 CHK(specint.quadrature_cdf[specint.quadrature_length-1] == 1.0); 274 } 275 276 /* Per SW spectral interval data */ 277 CHK(htgop_get_sw_spectral_intervals_wave_numbers(NULL, NULL) == RES_BAD_ARG); 278 CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, NULL) == RES_BAD_ARG); 279 CHK(htgop_get_sw_spectral_intervals_wave_numbers(NULL, &wnums) == RES_BAD_ARG); 280 CHK(htgop_get_sw_spectral_intervals_wave_numbers(htgop, &wnums) == RES_OK); 281 CHK(htgop_get_sw_spectral_interval(NULL, nspecints_sw, NULL) == RES_BAD_ARG); 282 CHK(htgop_get_sw_spectral_interval(htgop, nspecints_sw, NULL) == RES_BAD_ARG); 283 CHK(htgop_get_sw_spectral_interval(NULL, 0, NULL) == RES_BAD_ARG); 284 CHK(htgop_get_sw_spectral_interval(htgop, 0, NULL) == RES_BAD_ARG); 285 CHK(htgop_get_sw_spectral_interval(NULL, nspecints_sw, &specint) == RES_BAD_ARG); 286 CHK(htgop_get_sw_spectral_interval(htgop, nspecints_sw, &specint) == RES_BAD_ARG); 287 CHK(htgop_get_sw_spectral_interval(NULL, 0, &specint) == RES_BAD_ARG); 288 FOR_EACH(ispecint, 0, nspecints_sw) { 289 int i; 290 CHK(htgop_get_sw_spectral_interval(htgop, ispecint, &specint) == RES_OK); 291 CHK(specint.wave_numbers[0] == wnums[ispecint+0]); 292 CHK(specint.wave_numbers[1] == wnums[ispecint+1]); 293 294 /* Check the finding of a SW spectral interval */ 295 FOR_EACH(i, 0, 10) { 296 double wnum; 297 size_t id; 298 wnum = wnums[ispecint+0] 299 + rand_canonic() * (wnums[ispecint+1] - wnums[ispecint+0]); 300 CHK(htgop_find_sw_spectral_interval_id(htgop, wnum, &id) == RES_OK); 301 CHK(id == ispecint); 302 } 303 304 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 305 CHK(specint.wave_numbers[0] == dbl); 306 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 307 CHK(specint.wave_numbers[1] == dbl); 308 CHK(cstr_to_ulong(read_line(&rdr), &ul) == RES_OK); 309 CHK(specint.quadrature_length == ul); 310 311 FOR_EACH(iquad, 0, specint.quadrature_length) { 312 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 313 CHK(specint.quadrature_pdf[iquad] == dbl); 314 if(iquad == 0) { 315 CHK(specint.quadrature_cdf[iquad] == dbl); 316 } else { 317 const double pdf = 318 specint.quadrature_cdf[iquad] - specint.quadrature_cdf[iquad-1]; 319 CHK(eq_eps(pdf, dbl, 1.e-6)); 320 } 321 CHK(specint.quadrature_cdf[specint.quadrature_length-1] == 1.0); 322 } 323 } 324 325 /* Per layer LW ka data */ 326 FOR_EACH(ispecint, 0, nspecints_lw) { 327 CHK(htgop_get_lw_spectral_interval(htgop, ispecint, &specint) == RES_OK); 328 FOR_EACH(iquad, 0, specint.quadrature_length) { 329 FOR_EACH(ilay, 0, nlays) { 330 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 331 CHK(lay.lw_spectral_intervals_count == nspecints_lw); 332 CHK(lay.sw_spectral_intervals_count == nspecints_sw); 333 CHK(htgop_layer_get_lw_spectral_interval(&lay, ispecint, &lw_specint) == RES_OK); 334 CHK(lw_specint.quadrature_length == specint.quadrature_length); 335 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 336 CHK(lw_specint.ka_nominal[iquad] == dbl); 337 } 338 FOR_EACH(itab, 0, tab_len) { 339 FOR_EACH(ilay, 0, nlays) { 340 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 341 CHK(lay.lw_spectral_intervals_count == nspecints_lw); 342 CHK(lay.sw_spectral_intervals_count == nspecints_sw); 343 CHK(htgop_layer_get_lw_spectral_interval(&lay, ispecint, &lw_specint) == RES_OK); 344 CHK(htgop_layer_lw_spectral_interval_get_tab(&lw_specint, iquad, &lw_tab) == RES_OK); 345 CHK(lw_tab.tab_length == tab_len); 346 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 347 CHK(lw_tab.ka_tab[itab] == dbl); 348 CHK(lw_tab.x_h2o_tab[itab] == lay.x_h2o_tab[itab]); 349 } 350 } 351 } 352 } 353 354 /* Per layer SW ka data */ 355 FOR_EACH(ispecint, 0, nspecints_sw) { 356 CHK(htgop_get_sw_spectral_interval(htgop, ispecint, &specint) == RES_OK); 357 FOR_EACH(iquad, 0, specint.quadrature_length) { 358 FOR_EACH(ilay, 0, nlays) { 359 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 360 CHK(htgop_layer_get_sw_spectral_interval(&lay, ispecint, &sw_specint) == RES_OK); 361 CHK(sw_specint.quadrature_length == specint.quadrature_length); 362 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 363 CHK(sw_specint.ka_nominal[iquad] == dbl); 364 } 365 FOR_EACH(itab, 0, tab_len) { 366 FOR_EACH(ilay, 0, nlays) { 367 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 368 CHK(htgop_layer_get_sw_spectral_interval(&lay, ispecint, &sw_specint) == RES_OK); 369 CHK(htgop_layer_sw_spectral_interval_get_tab(&sw_specint, iquad, &sw_tab) == RES_OK); 370 CHK(sw_tab.tab_length == tab_len); 371 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 372 CHK(sw_tab.ka_tab[itab] == dbl); 373 CHK(sw_tab.x_h2o_tab[itab] == lay.x_h2o_tab[itab]); 374 } 375 } 376 } 377 } 378 379 /* Per layer SW ks data */ 380 FOR_EACH(ispecint, 0, nspecints_sw) { 381 CHK(htgop_get_sw_spectral_interval(htgop, ispecint, &specint) == RES_OK); 382 FOR_EACH(iquad, 0, specint.quadrature_length) { 383 FOR_EACH(ilay, 0, nlays) { 384 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 385 CHK(htgop_layer_get_sw_spectral_interval(&lay, ispecint, &sw_specint) == RES_OK); 386 CHK(sw_specint.quadrature_length == specint.quadrature_length); 387 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 388 CHK(sw_specint.ks_nominal[iquad] == dbl); 389 } 390 FOR_EACH(itab, 0, tab_len) { 391 FOR_EACH(ilay, 0, nlays) { 392 CHK(htgop_get_layer(htgop, ilay, &lay) == RES_OK); 393 CHK(htgop_layer_get_sw_spectral_interval(&lay, ispecint, &sw_specint) == RES_OK); 394 CHK(htgop_layer_sw_spectral_interval_get_tab(&sw_specint, iquad, &sw_tab) == RES_OK); 395 CHK(sw_tab.tab_length == tab_len); 396 CHK(cstr_to_double(read_line(&rdr), &dbl) == RES_OK); 397 CHK(sw_tab.ks_tab[itab] == dbl); 398 CHK(sw_tab.x_h2o_tab[itab] == lay.x_h2o_tab[itab]); 399 } 400 } 401 } 402 } 403 404 check_position_to_layer_id(htgop); 405 check_sw_specints(htgop, 1.e7/780.0, 1.e7/380.0); 406 check_lw_specints(htgop, 1.e7/100000, 1.e7/1000); 407 check_lw_specints(htgop, 1.e7/5600, 1.e7/3000); 408 check_lw_specints(htgop, 1.e7/9500, 1.e7/9300); 409 check_lw_specints(htgop, 1.e7/5000, 1.e7/5000); 410 411 CHK(fclose(fp) == 0); 412 reader_release(&rdr); 413 CHK(htgop_ref_put(htgop) == RES_OK); 414 check_memory_allocator(&allocator); 415 mem_shutdown_proxy_allocator(&allocator); 416 CHK(mem_allocated_size() == 0); 417 return 0; 418 } 419