1 /*
2 This file is part of Code_Saturne, a general-purpose CFD tool.
3
4 Copyright (C) 1998-2021 EDF S.A.
5
6 This program is free software; you can redistribute it and/or modify it under
7 the terms of the GNU General Public License as published by the Free Software
8 Foundation; either version 2 of the License, or (at your option) any later
9 version.
10
11 This program is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 details.
15
16 You should have received a copy of the GNU General Public License along with
17 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
18 Street, Fifth Floor, Boston, MA 02110-1301, USA. */
19
20 /*----------------------------------------------------------------------------*/
21
22 #include "cs_defs.h"
23
24 /*----------------------------------------------------------------------------
25 * Standard C library headers
26 *----------------------------------------------------------------------------*/
27
28 #include <assert.h>
29 #include <errno.h>
30 #include <ctype.h>
31 #include <stdio.h>
32 #include <stdarg.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <math.h>
36 #include <float.h>
37
38 #if defined(HAVE_MPI)
39 #include <mpi.h>
40 #endif
41
42 #if defined(HAVE_DLOPEN)
43 #include <dlfcn.h>
44 #endif
45
46 /*----------------------------------------------------------------------------
47 * Local headers
48 *----------------------------------------------------------------------------*/
49
50 #include "bft_mem.h"
51 #include "bft_error.h"
52 #include "bft_printf.h"
53
54 #include "cs_atmo_profile_std.h"
55 #include "cs_base.h"
56 #include "cs_boundary_conditions.h"
57 #include "cs_boundary_zone.h"
58 #include "cs_domain.h"
59 #include "cs_field.h"
60 #include "cs_field_default.h"
61 #include "cs_field_pointer.h"
62 #include "cs_halo.h"
63 #include "cs_log.h"
64 #include "cs_math.h"
65 #include "cs_mesh.h"
66 #include "cs_mesh_location.h"
67 #include "cs_mesh_quantities.h"
68 #include "cs_parall.h"
69 #include "cs_equation_iterative_solve.h"
70 #include "cs_physical_constants.h"
71 #include "cs_physical_model.h"
72 #include "cs_prototypes.h"
73 #include "cs_thermal_model.h"
74 #include "cs_turbulence_model.h"
75 #include "cs_volume_zone.h"
76 #include "cs_balance.h"
77 #include "cs_blas.h"
78 #include "cs_convection_diffusion.h"
79 #include "cs_parameters.h"
80 #include "cs_porous_model.h"
81 #include "cs_timer.h"
82 #include "cs_matrix_building.h"
83 #include "cs_sles.h"
84 #include "cs_sles_default.h"
85 #include "cs_face_viscosity.h"
86 #include "cs_divergence.h"
87 #include "cs_velocity_pressure.h"
88
89 /*----------------------------------------------------------------------------
90 * Header for the current file
91 *----------------------------------------------------------------------------*/
92
93 #include "cs_atmo.h"
94 #include "cs_atmo_aerosol.h"
95
96 /*----------------------------------------------------------------------------*/
97
98 BEGIN_C_DECLS
99
100 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
101
102 /*=============================================================================
103 * Local Macro Definitions
104 *============================================================================*/
105
106 /*=============================================================================
107 * Local Type Definitions
108 *============================================================================*/
109
110 /* global atmo options structure */
111 static cs_atmo_option_t _atmo_option = {
112 .syear = -1,
113 .squant = -1,
114 .shour = -1,
115 .smin = -1,
116 .ssec = -1.,
117 .longitude = 1e12, // TODO use cs_math_big_r
118 .latitude = 1e12,
119 .x_l93= 1e12,
120 .y_l93 = 1e12,
121 .domain_orientation = 0.,
122 .compute_z_ground = false,
123 .open_bcs_treatment = 0,
124 .sedimentation_model = 0,
125 .deposition_model = 0,
126 .nucleation_model = 0,
127 .subgrid_model = 0,
128 .distribution_model = 1, /* all or nothing */
129 .meteo_profile = 0, /* no meteo profile */
130 .meteo_file_name = NULL,
131 .meteo_dlmo = 0.,
132 .meteo_z0 = -1.,
133 .meteo_zref = -1.,
134 .meteo_zi = -1.,
135 .meteo_zu1 = -1.,
136 .meteo_zu2 = -1.,
137 .meteo_zt1 = -1.,
138 .meteo_zt2 = -1.,
139 .meteo_uref = -1.,
140 .meteo_u1 = -1.,
141 .meteo_u2 = -1.,
142 .meteo_ustar0 = -1.,
143 .meteo_wstar0 = -1.,
144 .meteo_angle = -1.,
145 .meteo_t0 = 284.15, /* 11 deg C */
146 .meteo_t1 = 0.,
147 .meteo_t2 = 0.,
148 .meteo_tstar = 0.,
149 .meteo_psea = 101325.,
150 .nbmetd = 0,
151 .nbmett = 0,
152 .nbmetm = 0,
153 .nbmaxt = 0,
154 .z_temp_met = NULL,
155 .time_met = NULL,
156 .hyd_p_met = NULL
157 };
158
159 /* global atmo constants structure */
160 static cs_atmo_constants_t _atmo_constants = {
161 .ps = 1.e5
162 };
163
164 /* atmo chemistry options structure */
165 static cs_atmo_chemistry_t _atmo_chem = {
166 .model = 0,
167 .n_species = 0,
168 .n_reactions = 0,
169 .chemistry_with_photolysis = true,
170 .aerosol_model = CS_ATMO_AEROSOL_OFF,
171 .frozen_gas_chem = false,
172 .init_gas_with_lib = false,
173 .init_aero_with_lib = false,
174 .n_layer = 0,
175 .n_size = 0,
176 .spack_file_name = NULL,
177 .species_to_scalar_id = NULL,
178 .species_to_field_id = NULL,
179 .molar_mass = NULL,
180 .chempoint = NULL,
181 .aero_file_name = NULL,
182 .chem_conc_file_name = NULL,
183 .aero_conc_file_name = NULL
184 };
185
186 /*============================================================================
187 * Static global variables
188 *============================================================================*/
189
190 cs_atmo_option_t *cs_glob_atmo_option = &_atmo_option;
191
192 cs_atmo_constants_t *cs_glob_atmo_constants = &_atmo_constants;
193
194 cs_atmo_chemistry_t *cs_glob_atmo_chemistry = &_atmo_chem;
195
196 static const char *cs_atmo_aerosol_type_enum_name[]
197 = {"CS_ATMO_AEROSOL_OFF",
198 "CS_ATMO_AEROSOL_SSH"};
199
200 static const char *cs_atmo_aerosol_type_name[]
201 = {N_("No atmospheric aerosol"),
202 N_("Atmospheric aerosol using external code SSH-aerosol")};
203
204 /*============================================================================
205 * Prototypes for functions intended for use only by Fortran wrappers.
206 * (descriptions follow, with function bodies).
207 *============================================================================*/
208
209 void
210 cs_f_atmo_get_meteo_file_name(int name_max,
211 const char **name,
212 int *name_len);
213
214 void
215 cs_f_atmo_get_chem_conc_file_name(int name_max,
216 const char **name,
217 int *name_len);
218
219 void
220 cs_f_atmo_get_aero_conc_file_name(int name_max,
221 const char **name,
222 int *name_len);
223
224 void
225 cs_f_atmo_get_pointers(cs_real_t **ps,
226 int **syear,
227 int **squant,
228 int **shour,
229 int **smin,
230 cs_real_t **ssec,
231 cs_real_t **longitude,
232 cs_real_t **latitude,
233 cs_real_t **x_l93,
234 cs_real_t **y_l93,
235 bool **compute_z_ground,
236 int **open_bcs_treatment,
237 int **sedimentation_model,
238 int **deposition_model,
239 int **nucleation_model,
240 int **subgrid_model,
241 int **distribution_model,
242 int **model,
243 int **n_species,
244 int **n_reactions,
245 bool **chemistry_with_photolysis,
246 cs_atmo_aerosol_type_t **aerosol_model,
247 bool **frozen_gas_chem,
248 bool **init_gas_with_lib,
249 bool **init_aero_with_lib,
250 int **n_layer,
251 int **n_size,
252 int **meteo_profile,
253 int **nbmetd,
254 int **nbmett,
255 int **nbmetm,
256 int **nbmaxt,
257 cs_real_t **meteo_zi);
258
259 void
260 cs_f_atmo_arrays_get_pointers(cs_real_t **z_temp_met,
261 cs_real_t **time_met,
262 cs_real_t **hyd_p_met,
263 int dim_hyd_p_met[2]);
264
265 void
266 cs_f_atmo_chem_arrays_get_pointers(int **species_to_scalar_id,
267 cs_real_t **molar_mass,
268 int **chempoint);
269
270 void
271 cs_f_atmo_chem_initialize_species_to_fid(int *species_fid);
272
273 void
274 cs_f_atmo_chem_finalize(void);
275
276 /*============================================================================
277 * Fortran wrapper function definitions
278 *============================================================================*/
279
280 /*----------------------------------------------------------------------------
281 * Return the name of the meteo file
282 *
283 * This function is intended for use by Fortran wrappers.
284 *
285 * parameters:
286 * name_max <-- maximum name length
287 * name --> pointer to associated length
288 * name_len --> length of associated length
289 *----------------------------------------------------------------------------*/
290
291 void
cs_f_atmo_get_meteo_file_name(int name_max,const char ** name,int * name_len)292 cs_f_atmo_get_meteo_file_name(int name_max,
293 const char **name,
294 int *name_len)
295 {
296 *name = _atmo_option.meteo_file_name;
297 *name_len = strlen(*name);
298
299 if (*name_len > name_max) {
300 bft_error
301 (__FILE__, __LINE__, 0,
302 _("Error retrieving meteo file (\"%s\"):\n"
303 "Fortran caller name length (%d) is too small for name \"%s\"\n"
304 "(of length %d)."),
305 _atmo_option.meteo_file_name, name_max, *name, *name_len);
306 }
307 }
308
309 /*----------------------------------------------------------------------------
310 * Return the name of the chemistry concentration file
311 *
312 * This function is intended for use by Fortran wrappers.
313 *
314 * parameters:
315 * name_max <-- maximum name length
316 * name --> pointer to associated length
317 * name_len --> length of associated length
318 *----------------------------------------------------------------------------*/
319
320 void
cs_f_atmo_get_chem_conc_file_name(int name_max,const char ** name,int * name_len)321 cs_f_atmo_get_chem_conc_file_name(int name_max,
322 const char **name,
323 int *name_len)
324 {
325 *name = _atmo_chem.chem_conc_file_name;
326 *name_len = strlen(*name);
327
328 if (*name_len > name_max) {
329 bft_error
330 (__FILE__, __LINE__, 0,
331 _("Error retrieving chemistry concentration file (\"%s\"):\n"
332 "Fortran caller name length (%d) is too small for name \"%s\"\n"
333 "(of length %d)."),
334 _atmo_chem.chem_conc_file_name, name_max, *name, *name_len);
335 }
336 }
337
338 /*----------------------------------------------------------------------------
339 * Return the name of the aerosol concentration file
340 *
341 * This function is intended for use by Fortran wrappers.
342 *
343 * parameters:
344 * name_max <-- maximum name length
345 * name --> pointer to associated length
346 * name_len --> length of associated length
347 *----------------------------------------------------------------------------*/
348
349 void
cs_f_atmo_get_aero_conc_file_name(int name_max,const char ** name,int * name_len)350 cs_f_atmo_get_aero_conc_file_name(int name_max,
351 const char **name,
352 int *name_len)
353 {
354 *name = _atmo_chem.aero_conc_file_name;
355 *name_len = strlen(*name);
356
357 if (*name_len > name_max) {
358 bft_error
359 (__FILE__, __LINE__, 0,
360 _("Error retrieving chemistry concentration file (\"%s\"):\n"
361 "Fortran caller name length (%d) is too small for name \"%s\"\n"
362 "(of length %d)."),
363 _atmo_chem.aero_conc_file_name, name_max, *name, *name_len);
364 }
365 }
366
367 /*----------------------------------------------------------------------------
368 * Access pointers for Fortran mapping.
369 *----------------------------------------------------------------------------*/
370
371 void
cs_f_atmo_get_pointers(cs_real_t ** ps,int ** syear,int ** squant,int ** shour,int ** smin,cs_real_t ** ssec,cs_real_t ** longitude,cs_real_t ** latitude,cs_real_t ** x_l93,cs_real_t ** y_l93,bool ** compute_z_ground,int ** open_bcs_treatment,int ** sedimentation_model,int ** deposition_model,int ** nucleation_model,int ** subgrid_model,int ** distribution_model,int ** model,int ** n_species,int ** n_reactions,bool ** chemistry_with_photolysis,cs_atmo_aerosol_type_t ** aerosol_model,bool ** frozen_gas_chem,bool ** init_gas_with_lib,bool ** init_aero_with_lib,int ** n_layer,int ** n_size,int ** meteo_profile,int ** nbmetd,int ** nbmett,int ** nbmetm,int ** nbmaxt,cs_real_t ** meteo_zi)372 cs_f_atmo_get_pointers(cs_real_t **ps,
373 int **syear,
374 int **squant,
375 int **shour,
376 int **smin,
377 cs_real_t **ssec,
378 cs_real_t **longitude,
379 cs_real_t **latitude,
380 cs_real_t **x_l93,
381 cs_real_t **y_l93,
382 bool **compute_z_ground,
383 int **open_bcs_treatment,
384 int **sedimentation_model,
385 int **deposition_model,
386 int **nucleation_model,
387 int **subgrid_model,
388 int **distribution_model,
389 int **model,
390 int **n_species,
391 int **n_reactions,
392 bool **chemistry_with_photolysis,
393 cs_atmo_aerosol_type_t **aerosol_model,
394 bool **frozen_gas_chem,
395 bool **init_gas_with_lib,
396 bool **init_aero_with_lib,
397 int **n_layer,
398 int **n_size,
399 int **meteo_profile,
400 int **nbmetd,
401 int **nbmett,
402 int **nbmetm,
403 int **nbmaxt,
404 cs_real_t **meteo_zi)
405 {
406 *ps = &(_atmo_constants.ps);
407 *syear = &(_atmo_option.syear);
408 *squant = &(_atmo_option.squant);
409 *shour = &(_atmo_option.shour);
410 *smin = &(_atmo_option.smin);
411 *ssec = &(_atmo_option.ssec);
412 *longitude = &(_atmo_option.longitude);
413 *latitude = &(_atmo_option.latitude);
414 *x_l93 = &(_atmo_option.x_l93);
415 *y_l93 = &(_atmo_option.y_l93);
416 *compute_z_ground = &(_atmo_option.compute_z_ground);
417 *open_bcs_treatment = &(_atmo_option.open_bcs_treatment);
418 *sedimentation_model = &(_atmo_option.sedimentation_model);
419 *deposition_model = &(_atmo_option.deposition_model);
420 *nucleation_model = &(_atmo_option.nucleation_model);
421 *subgrid_model = &(_atmo_option.subgrid_model);
422 *distribution_model = &(_atmo_option.distribution_model);
423 *meteo_profile = &(_atmo_option.meteo_profile);
424 *nbmetd = &(_atmo_option.nbmetd);
425 *nbmett = &(_atmo_option.nbmett);
426 *nbmetm = &(_atmo_option.nbmetm);
427 *nbmaxt = &(_atmo_option.nbmaxt);
428 *meteo_zi = &(_atmo_option.meteo_zi);
429 *model = &(_atmo_chem.model);
430 *n_species = &(_atmo_chem.n_species);
431 *n_reactions = &(_atmo_chem.n_reactions);
432 *chemistry_with_photolysis = &(_atmo_chem.chemistry_with_photolysis);
433 *aerosol_model = &(_atmo_chem.aerosol_model);
434 *frozen_gas_chem = &(_atmo_chem.frozen_gas_chem);
435 *init_gas_with_lib = &(_atmo_chem.init_gas_with_lib);
436 *init_aero_with_lib = &(_atmo_chem.init_aero_with_lib);
437 *n_layer = &(_atmo_chem.n_layer);
438 *n_size = &(_atmo_chem.n_size);
439 }
440
441 void
cs_f_atmo_chem_arrays_get_pointers(int ** species_to_scalar_id,cs_real_t ** molar_mass,int ** chempoint)442 cs_f_atmo_chem_arrays_get_pointers(int **species_to_scalar_id,
443 cs_real_t **molar_mass,
444 int **chempoint)
445 {
446 if (_atmo_chem.species_to_scalar_id == NULL)
447 BFT_MALLOC(_atmo_chem.species_to_scalar_id, _atmo_chem.n_species, int);
448 if (_atmo_chem.species_to_field_id == NULL)
449 BFT_MALLOC(_atmo_chem.species_to_field_id, _atmo_chem.n_species, int);
450 if (_atmo_chem.molar_mass == NULL)
451 BFT_MALLOC(_atmo_chem.molar_mass, _atmo_chem.n_species, cs_real_t);
452 if (_atmo_chem.chempoint == NULL)
453 BFT_MALLOC(_atmo_chem.chempoint, _atmo_chem.n_species, int);
454
455 *species_to_scalar_id = (_atmo_chem.species_to_scalar_id);
456 *molar_mass = (_atmo_chem.molar_mass);
457 *chempoint = (_atmo_chem.chempoint);
458 }
459
460 void
cs_f_atmo_arrays_get_pointers(cs_real_t ** z_temp_met,cs_real_t ** time_met,cs_real_t ** hyd_p_met,int dim_hyd_p_met[2])461 cs_f_atmo_arrays_get_pointers(cs_real_t **z_temp_met,
462 cs_real_t **time_met,
463 cs_real_t **hyd_p_met,
464 int dim_hyd_p_met[2])
465 {
466 if (_atmo_option.z_temp_met == NULL)
467 BFT_MALLOC(_atmo_option.z_temp_met, _atmo_option.nbmaxt, cs_real_t);
468 if (_atmo_option.time_met == NULL)
469 BFT_MALLOC(_atmo_option.time_met, _atmo_option.nbmetm, cs_real_t);
470 if (_atmo_option.hyd_p_met == NULL)
471 BFT_MALLOC(_atmo_option.hyd_p_met,
472 _atmo_option.nbmetm*_atmo_option.nbmaxt, cs_real_t);
473
474 *hyd_p_met = _atmo_option.hyd_p_met;
475 dim_hyd_p_met[0] = _atmo_option.nbmaxt;
476 dim_hyd_p_met[1] = _atmo_option.nbmetm;
477
478 *z_temp_met = _atmo_option.z_temp_met;
479 *time_met = _atmo_option.time_met;
480 }
481
482 void
cs_f_atmo_chem_initialize_species_to_fid(int * species_fid)483 cs_f_atmo_chem_initialize_species_to_fid(int *species_fid)
484 {
485 assert(species_fid != NULL);
486 assert(_atmo_chem.species_to_field_id != NULL);
487
488 for (int i = 0; i < _atmo_chem.n_species; i++)
489 _atmo_chem.species_to_field_id[i] = species_fid[i];
490 }
491
492 void
cs_f_atmo_chem_finalize(void)493 cs_f_atmo_chem_finalize(void)
494 {
495 if (_atmo_chem.aerosol_model != CS_ATMO_AEROSOL_OFF)
496 cs_atmo_aerosol_finalize();
497
498 BFT_FREE(_atmo_chem.species_to_scalar_id);
499 BFT_FREE(_atmo_chem.species_to_field_id);
500 BFT_FREE(_atmo_chem.molar_mass);
501 BFT_FREE(_atmo_chem.chempoint);
502 BFT_FREE(_atmo_chem.spack_file_name);
503 BFT_FREE(_atmo_chem.aero_file_name);
504 BFT_FREE(_atmo_chem.chem_conc_file_name);
505 }
506
507 /*============================================================================
508 * Private function definitions
509 *============================================================================*/
510
511 /*----------------------------------------------------------------------------
512 * Convert string to lower case
513 *----------------------------------------------------------------------------*/
514
515 static void
_strtolower(char * dest,const char * src)516 _strtolower(char *dest,
517 const char *src)
518 {
519 char *_dest = dest;
520 while (*src) {
521 *_dest = tolower(*src);
522 src++;
523 _dest++;
524 }
525 }
526
527 /*----------------------------------------------------------------------------*/
528 /*!
529 * \brief This auxiliary function solve a Poisson equation for hydrostatic
530 * pressure :
531 * \f$ \divs ( \grad P ) = \divs ( f ) \f$
532 * \param[in] f_ext external forcing
533 * \param[in] dfext external forcing increment
534 * \param[in] name field name
535 * \param[in] min_z Minimum altitude of the domain
536 * \param[in] p_ground Pressure at the minimum altitude
537 * \param[in,out] i_massflux Internal mass flux
538 * \param[in,out] b_massflux Boundary mass flux
539 * \param[in,out] i_viscm Internal face viscosity
540 * \param[in,out] i_viscm Boundary face viscosity
541 * \param[in,out] dam Working array
542 * \param[in,out] xam Working array
543 * \param[in,out] dpvar Pressure increment
544 * \param[in,out] rhs Working array
545 */
546 /*----------------------------------------------------------------------------*/
547
548 static void
_hydrostatic_pressure_compute(cs_real_3_t f_ext[],cs_real_3_t dfext[],cs_real_t pvar[],const char * name,cs_real_t min_z,cs_real_t p_ground,cs_real_t i_massflux[],cs_real_t b_massflux[],cs_real_t i_viscm[],cs_real_t b_viscm[],cs_real_t dam[],cs_real_t xam[],cs_real_t dpvar[],cs_real_t rhs[])549 _hydrostatic_pressure_compute(cs_real_3_t f_ext[],
550 cs_real_3_t dfext[],
551 cs_real_t pvar[],
552 const char *name,
553 cs_real_t min_z,
554 cs_real_t p_ground,
555 cs_real_t i_massflux[],
556 cs_real_t b_massflux[],
557 cs_real_t i_viscm[],
558 cs_real_t b_viscm[],
559 cs_real_t dam[],
560 cs_real_t xam[],
561 cs_real_t dpvar[],
562 cs_real_t rhs[])
563
564 {
565 /* Local variables */
566 cs_domain_t *domain = cs_glob_domain;
567 cs_mesh_t *m = domain->mesh;
568 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
569 const cs_real_3_t *restrict b_face_cog
570 = (const cs_real_3_t *restrict)mq->b_face_cog;
571 cs_field_t *f = cs_field_by_name(name);
572 int f_id = f->id;
573 int niterf;
574
575 cs_var_cal_opt_t vcopt;
576 cs_field_get_key_struct(f, cs_field_key_id("var_cal_opt"), &vcopt);
577
578 /*==========================================================================
579 * 0. Initialization
580 *==========================================================================*/
581
582 /* solving info */
583 cs_solving_info_t sinfo;
584 int key_sinfo_id = cs_field_key_id("solving_info");
585 if (f_id > -1) {
586 f = cs_field_by_id(f_id);
587 cs_field_get_key_struct(f, key_sinfo_id, &sinfo);
588 }
589
590 /* Symmetric matrix, except if advection */
591 int isym = 1;
592 bool symmetric = true;
593
594 /* Matrix block size */
595 cs_lnum_t eb_size[4], db_size[4];
596 int ibsize = 1, iesize = 1;
597
598 db_size[0] = ibsize;
599 db_size[1] = ibsize;
600 db_size[2] = ibsize;
601 db_size[3] = ibsize*ibsize;
602
603 eb_size[0] = iesize;
604 eb_size[1] = iesize;
605 eb_size[2] = iesize;
606 eb_size[3] = iesize*iesize;
607
608 cs_real_3_t *next_fext;
609 BFT_MALLOC(next_fext, m->n_cells_with_ghosts, cs_real_3_t);
610 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
611 next_fext[cell_id][0] = f_ext[cell_id][0] + dfext[cell_id][0];
612 next_fext[cell_id][1] = f_ext[cell_id][1] + dfext[cell_id][1];
613 next_fext[cell_id][2] = f_ext[cell_id][2] + dfext[cell_id][2];
614 }
615
616 /* --> Handle parallelism and periodicity */
617 if (cs_glob_rank_id >= 0 || cs_glob_mesh->n_init_perio > 0)
618 cs_mesh_sync_var_vect((cs_real_t *)next_fext);
619
620 /* Boundary conditions
621 *====================*/
622
623 vcopt.ndircl = 0;
624
625 /* To solve hydrostatic pressure:
626 * p_ground on the lowest face, homogeneous Neumann everywhere else */
627 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
628 cs_real_t hint = 1. / mq->b_dist[face_id];
629 cs_real_t qimp = 0.;
630 if ((b_face_cog[face_id][2] - min_z) < 0.005) {//FIXME dimensionless value
631 cs_real_t pimp = p_ground;
632 cs_boundary_conditions_set_dirichlet_scalar(&(f->bc_coeffs->a[face_id]),
633 &(f->bc_coeffs->af[face_id]),
634 &(f->bc_coeffs->b[face_id]),
635 &(f->bc_coeffs->bf[face_id]),
636 pimp,
637 hint,
638 cs_math_big_r);
639 vcopt.ndircl = 1;
640 }
641 else {
642 cs_boundary_conditions_set_neumann_scalar(&(f->bc_coeffs->a[face_id]),
643 &(f->bc_coeffs->af[face_id]),
644 &(f->bc_coeffs->b[face_id]),
645 &(f->bc_coeffs->bf[face_id]),
646 qimp,
647 hint);
648 }
649 }
650
651 cs_parall_max(1, CS_INT_TYPE, &(vcopt.ndircl));
652 cs_real_t *rovsdt;
653 BFT_MALLOC(rovsdt, m->n_cells_with_ghosts, cs_real_t);
654
655 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells_with_ghosts; cell_id++)
656 rovsdt[cell_id] = 0.;
657
658 /* Faces viscosity */
659 cs_real_t *c_visc;
660 BFT_MALLOC(c_visc, m->n_cells_with_ghosts, cs_real_t);
661
662 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells_with_ghosts; cell_id++)
663 c_visc[cell_id] = 1.;
664
665 cs_face_viscosity(m, mq, vcopt.imvisf, c_visc, i_viscm, b_viscm);
666
667 cs_matrix_wrapper_scalar(vcopt.iconv,
668 vcopt.idiff,
669 vcopt.ndircl,
670 isym,
671 vcopt.thetav,
672 0,
673 f->bc_coeffs->b,
674 f->bc_coeffs->bf,
675 rovsdt,
676 i_massflux,
677 b_massflux,
678 i_viscm,
679 b_viscm,
680 NULL,
681 dam,
682 xam);
683
684 /* Right hand side
685 *================*/
686
687 cs_ext_force_flux(m,
688 mq,
689 1, /* init */
690 vcopt.nswrgr,
691 next_fext,
692 f->bc_coeffs->bf,
693 i_massflux,
694 b_massflux,
695 i_viscm,
696 b_viscm,
697 c_visc,
698 c_visc,
699 c_visc);
700
701 cs_real_t *divergfext;
702 BFT_MALLOC(divergfext, m->n_cells_with_ghosts, cs_real_t);
703
704 cs_divergence(m,
705 1, /* init */
706 i_massflux,
707 b_massflux,
708 divergfext);
709
710 /* --- Right hand side residual */
711 cs_real_t rnorm = sqrt(cs_gdot(m->n_cells, divergfext, divergfext));
712 cs_real_t residu = rnorm;
713
714 /* Initial Right-Hand-Side */
715 cs_diffusion_potential(-1,//FIXME why not f_id?
716 m,
717 mq,
718 1, /* init */
719 1, /* inc */
720 vcopt.imrgra,
721 1, /* iccocg */
722 vcopt.nswrgr,
723 vcopt.imligr,
724 1, /* iphydp */
725 vcopt.iwgrec,
726 vcopt.iwarni,
727 vcopt.epsrgr,
728 vcopt.climgr,
729 next_fext,
730 pvar,
731 f->bc_coeffs->a,
732 f->bc_coeffs->b,
733 f->bc_coeffs->af,
734 f->bc_coeffs->bf,
735 i_viscm,
736 b_viscm,
737 c_visc,
738 rhs);
739
740 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++)
741 rhs[cell_id] = - divergfext[cell_id] - rhs[cell_id];
742
743 cs_real_t ressol = 0;
744 for (int sweep = 0;
745 sweep < vcopt.nswrsm && residu > (rnorm * vcopt.epsrsm);
746 sweep++) {
747
748 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++)
749 dpvar[cell_id] = 0.;
750
751 ressol = residu;
752 cs_sles_solve_native(f_id,
753 "",
754 symmetric,
755 db_size,
756 eb_size,
757 dam,
758 xam,
759 vcopt.epsilo,
760 rnorm,
761 &niterf,
762 &ressol,
763 rhs,
764 dpvar);
765
766 /* Update variable and right-hand-side */
767 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++)
768 pvar[cell_id] += dpvar[cell_id];
769
770
771 cs_diffusion_potential(-1,//FIXME why not f_id?
772 m,
773 mq,
774 1, /* init */
775 1, /* inc */
776 vcopt.imrgra,
777 1, /* iccocg */
778 vcopt.nswrgr,
779 vcopt.imligr,
780 1, /* iphydp */
781 vcopt.iwgrec,
782 vcopt.iwarni,
783 vcopt.epsrgr,
784 vcopt.climgr,
785 next_fext,
786 pvar,
787 f->bc_coeffs->a,
788 f->bc_coeffs->b,
789 f->bc_coeffs->af,
790 f->bc_coeffs->bf,
791 i_viscm,
792 b_viscm,
793 c_visc,
794 rhs);
795
796 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++)
797 rhs[cell_id] = - divergfext[cell_id] - rhs[cell_id];
798
799 /* --- Convergence test */
800 residu = sqrt(cs_gdot(m->n_cells, rhs, rhs));
801
802 /* Writing */
803 if (vcopt.iwarni >= 2) {
804 bft_printf("%s: CV_DIF_TS, IT: %d, Res: %12.5e, Norm: %12.5e\n",
805 name, sweep, residu, rnorm);
806 bft_printf("%s: Current reconstruction sweep: %d, "
807 "Iterations for solver: %d\n", name, sweep, niterf);
808 }
809
810 }
811
812 cs_sles_free_native(f_id, "");
813
814 BFT_FREE(divergfext);
815 BFT_FREE(next_fext);
816 BFT_FREE(rovsdt);
817 BFT_FREE(c_visc);
818 }
819
820 /*----------------------------------------------------------------------------*/
821 /*!
822 * \brief Convert coordinates from Lambert-93 to WGS84 inside atmo options
823 */
824 /*----------------------------------------------------------------------------*/
825
826 static void
_convert_from_l93_to_wgs84(void)827 _convert_from_l93_to_wgs84(void)
828 {
829 //computation from https://georezo.net/forum/viewtopic.php?id=94465
830 cs_real_t c = 11754255.426096; // projection constant
831 cs_real_t e = 0.0818191910428158; // ellipsoid excentricity
832 cs_real_t n = 0.725607765053267; // projection exponent
833 cs_real_t xs = 700000; // projection's pole x-coordinate
834 cs_real_t ys = 12655612.049876; // projection's pole y-coordinate
835 cs_real_t a = (log(c/(sqrt( cs_math_pow2(cs_glob_atmo_option->x_l93-xs)
836 + cs_math_pow2(cs_glob_atmo_option->y_l93-ys))))/n);
837
838 double t1 = a + e*atanh(e*(tanh(a+e*atanh(e*sin(1)))));
839 double t2 = e*tanh(a+e*atanh(e*(tanh(t1))));
840 double t3 = a+e*atanh(e*tanh(a+e*atanh(e*tanh(a+e*atanh(t2)))));
841
842 cs_glob_atmo_option->longitude = ((atan(-(cs_glob_atmo_option->x_l93-xs)
843 /(cs_glob_atmo_option->y_l93-ys)))/n
844 + 3*cs_math_pi/180)/cs_math_pi*180;
845 cs_glob_atmo_option->latitude
846 = asin(tanh((log(c/sqrt( cs_math_pow2(cs_glob_atmo_option->x_l93-xs)
847 + cs_math_pow2(cs_glob_atmo_option->y_l93-ys)))/n)
848 +e*atanh(e*tanh(t3))))/cs_math_pi*180;
849 }
850
851 /*----------------------------------------------------------------------------*/
852 /*!
853 * \brief Convert coordinates from WGS84 to Lambert-93 inside atmo options
854 */
855 /*----------------------------------------------------------------------------*/
856
857 static void
_convert_from_wgs84_to_l93(void)858 _convert_from_wgs84_to_l93(void)
859 {
860 //computation from https://georezo.net/forum/viewtopic.php?id=94465
861 cs_real_t c = 11754255.426096; // projection constant
862 cs_real_t e = 0.0818191910428158; // ellipsoid excentricity
863 cs_real_t n = 0.725607765053267; // projection exponent
864 cs_real_t xs = 700000; // projection's pole x-coordinate
865 cs_real_t ys = 12655612.049876; // projection's pole y-coordinate
866
867 cs_real_t lat_rad= cs_glob_atmo_option->latitude*cs_math_pi/180; //latitude in rad
868 cs_real_t lat_iso= atanh(sin(lat_rad))-e*atanh(e*sin(lat_rad)); // isometric latitude
869
870 cs_glob_atmo_option->x_l93= ((c*exp(-n*(lat_iso)))
871 *sin(n*(cs_glob_atmo_option->longitude-3)
872 *cs_math_pi/180)+xs);
873 cs_glob_atmo_option->y_l93= (ys-(c*exp(-n*(lat_iso)))
874 *cos(n*(cs_glob_atmo_option->longitude-3)
875 *cs_math_pi/180));
876 }
877
878 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
879
880 /*============================================================================
881 * Public function definitions
882 *============================================================================*/
883
884 /*----------------------------------------------------------------------------*/
885 /*!
886 * \brief Initialize meteo profiles if no meteo file is given.
887 */
888 /*----------------------------------------------------------------------------*/
889
890 void
cs_atmo_init_meteo_profiles(void)891 cs_atmo_init_meteo_profiles(void)
892 {
893 /* Some turbulence constants */
894 cs_real_t kappa = cs_turb_xkappa;
895
896 cs_atmo_option_t *aopt = &_atmo_option;
897 cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();
898
899 /* potential temp at ref */
900 cs_real_t pref = cs_glob_atmo_constants->ps;
901 cs_real_t rair = phys_pro->r_pg_cnst;
902 cs_real_t cp0 = phys_pro->cp0;
903 cs_real_t rscp = rair/cp0;
904 cs_real_t g = cs_math_3_norm(cs_glob_physical_constants->gravity);
905 if (g <= 0.)
906 bft_error(__FILE__,__LINE__, 0,
907 _("Atmo meteo profiles: gravity must not be 0.\n"));
908
909 cs_real_t theta0 = aopt->meteo_t0 * pow(pref/ aopt->meteo_psea, rscp);
910
911 /* Reference fluid properties set from meteo values */
912 phys_pro->p0 = aopt->meteo_psea;
913 phys_pro->t0 = aopt->meteo_t0; /* ref temp T0 */
914 phys_pro->ro0 = phys_pro->p0/(rair * aopt->meteo_t0); /* ref density T0 */
915
916 cs_real_t z0 = aopt->meteo_z0;
917 cs_real_t zref = aopt->meteo_zref;
918 if (aopt->meteo_ustar0 < 0. && aopt->meteo_uref < 0.)
919 bft_error(__FILE__,__LINE__, 0,
920 _("Atmo meteo profiles: meteo_ustar0 or meteo_uref.\n"));
921
922 /* Recompute LMO inverse */
923 if (aopt->meteo_ustar0 >= 0. && aopt->meteo_uref >= 0.) {
924
925 /* U+ */
926 cs_real_t up = aopt->meteo_uref / aopt->meteo_ustar0;
927
928 /* U+ if neutral, dlmo = 0 */
929 cs_real_t up_l = cs_mo_psim(zref + z0,
930 z0,
931 0.) / kappa;
932
933 cs_real_t dlmo = 0.;
934 cs_real_t error = up_l - up;
935
936 /* Dichotomy */
937 cs_real_t dl_min = -1.e6;
938 cs_real_t dl_max = 1.e6;
939 cs_real_t tol = 1e-6;
940 int it;
941 int it_max = 1000;
942 for (it = 0;
943 it < it_max && CS_ABS(error) > tol && 0.5*(dl_max - dl_min) > tol;
944 it++) {
945 cs_real_t dl_mid = 0.5 * (dl_min + dl_max);
946
947 cs_real_t error_min = cs_mo_psim(zref + z0,
948 z0,
949 dl_min) / kappa - up;
950 cs_real_t error_mid = cs_mo_psim(zref + z0,
951 z0,
952 dl_mid) / kappa - up;
953
954 /* The solution is between min and mid */
955 if (error_min * error_mid < 0) {
956 dl_max = dl_mid;
957 if (CS_ABS(error_min) < CS_ABS(error_mid)) {
958 dlmo = dl_min;
959 error = error_min;
960 }
961 else {
962 dlmo = dl_mid;
963 error = error_mid;
964 }
965 }
966 /* The solution is between mid and max */
967 else {
968 cs_real_t error_max = cs_mo_psim(zref + z0,
969 z0,
970 dl_max) / kappa - up;
971 dl_min = dl_mid;
972 if (CS_ABS(error_mid) < CS_ABS(error_max)) {
973 dlmo = dl_mid;
974 error = error_mid;
975 }
976 else {
977 dlmo = dl_max;
978 error = error_max;
979 }
980 }
981 #if 0
982 bft_printf("IT %d: dlmo = %f, error = %f\n",it, dlmo, error);
983 #endif
984 }
985
986 if (it == it_max)
987 bft_printf("Warning: meteo preprocessor did not converge to find inverse\n"
988 " of LMO length, current value is %f.\n"
989 "Please, check reference velocity, reference altitude and ustar\n",
990 dlmo);
991
992 aopt->meteo_dlmo = dlmo;
993 }
994
995 /* Compute ground friction velocity from dlmo and uref */
996 if (aopt->meteo_ustar0 < 0.)
997 aopt->meteo_ustar0 =
998 aopt->meteo_uref * kappa
999 / cs_mo_psim(zref + z0,
1000 z0,
1001 aopt->meteo_dlmo);
1002
1003 /* Compute uref from ground friction velocity and dlmo */
1004 if (aopt->meteo_uref < 0. && zref > 0.)
1005 aopt->meteo_uref =
1006 aopt->meteo_ustar0 / kappa
1007 * cs_mo_psim(zref + z0,
1008 z0,
1009 aopt->meteo_dlmo);
1010
1011 /* LMO inverse, ustar at ground */
1012 cs_real_t dlmo = aopt->meteo_dlmo;
1013 cs_real_t ustar0 = aopt->meteo_ustar0;
1014
1015 /* Friction temperature */
1016 aopt->meteo_tstar = cs_math_pow2(ustar0) * theta0 * dlmo / (kappa * g);
1017
1018 /* Center of the domain */
1019 /* if neither latitude/longitude nor lambert coordinates are given */
1020 if (((aopt->latitude > 0.5 * cs_math_big_r || aopt->longitude > 0.5 * cs_math_big_r)
1021 && (aopt->x_l93 > 0.5 * cs_math_big_r || aopt->y_l93 > 0.5 * cs_math_big_r))) {
1022 bft_printf("Neither latitude nor center in Lambert-93 was given \n");
1023 bft_printf("It is set to Paris values \n");
1024 aopt->latitude = 45.44;
1025 aopt->longitude = 4.39;
1026 _convert_from_wgs84_to_l93();
1027 }
1028 /* else if latitude/longitude are given */
1029 else if (aopt->x_l93 > 0.5 * cs_math_big_r || aopt->y_l93 > 0.5 * cs_math_big_r) {
1030 bft_printf("Latitude and longitude were given, Lambert center's coordinates"
1031 "are automatically computed\n");
1032 _convert_from_wgs84_to_l93();
1033 }
1034 /* All other cases */
1035 else{
1036 bft_printf("Lambert coordinates were given, latitude"
1037 "and longitude are automatically computed\n");
1038 _convert_from_l93_to_wgs84();
1039 }
1040
1041 /* BL height according to Marht 1982 formula */
1042 /* value of C=0.2, 0.185, 0.06, 0.14, 0.07, 0.04 */
1043 cs_real_t zi_coef = 0.2;
1044 cs_real_t corio_f = 4. * cs_math_pi / 86164.
1045 * sin(aopt->latitude * cs_math_pi / 180.);
1046 aopt->meteo_zi = zi_coef * ustar0 / fabs(corio_f);
1047
1048 /* Force the computation of z_ground */
1049 aopt->compute_z_ground = true;
1050 }
1051
1052 /*----------------------------------------------------------------------------*/
1053 /*!
1054 * \brief Compute meteo profiles if no meteo file is given.
1055 */
1056 /*----------------------------------------------------------------------------*/
1057
1058 void
cs_atmo_compute_meteo_profiles(void)1059 cs_atmo_compute_meteo_profiles(void)
1060 {
1061 cs_domain_t *domain = cs_glob_domain;
1062 cs_mesh_t *m = domain->mesh;
1063 cs_mesh_quantities_t *mq = domain->mesh_quantities;
1064
1065 const cs_real_3_t *restrict cell_cen
1066 = (const cs_real_3_t *restrict)mq->cell_cen;
1067
1068 /* In the log */
1069 bft_printf(" Computing meteo profiles from CS\n\n");
1070
1071 /* Get fields */
1072 cs_real_t *cpro_met_potemp = cs_field_by_name("meteo_pot_temperature")->val;
1073 cs_real_3_t *cpro_met_vel =
1074 (cs_real_3_t *) (cs_field_by_name("meteo_velocity")->val);
1075 cs_real_t *cpro_met_k = cs_field_by_name("meteo_tke")->val;
1076 cs_real_t *cpro_met_eps = cs_field_by_name("meteo_eps")->val;
1077
1078 /* Some turbulence constants */
1079 cs_real_t kappa = cs_turb_xkappa;
1080 cs_real_t cmu = cs_turb_cmu;
1081
1082 cs_atmo_option_t *aopt = &_atmo_option;
1083 cs_real_t z0 = aopt->meteo_z0;
1084
1085 const cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();
1086 /* potential temp at ref */
1087 cs_real_t pref = cs_glob_atmo_constants->ps;
1088 cs_real_t rair = phys_pro->r_pg_cnst;
1089 cs_real_t cp0 = phys_pro->cp0;
1090 cs_real_t rscp = rair/cp0;
1091 cs_real_t theta0 = aopt->meteo_t0 * pow(pref/ aopt->meteo_psea, rscp);
1092
1093 /* LMO inverse, ustar at ground */
1094 cs_real_t dlmo = aopt->meteo_dlmo;
1095 cs_real_t ustar0 = aopt->meteo_ustar0;
1096 cs_real_t angle = aopt->meteo_angle;
1097
1098 /* Friction temperature */
1099 cs_real_t tstar = aopt->meteo_tstar;
1100
1101 /* Variables used for clipping */
1102 cs_real_t ri_max = cs_math_big_r;
1103 cs_real_t *dlmo_var = NULL;
1104 cs_real_t z_lim = cs_math_big_r;
1105 cs_real_t u_met_min= cs_math_big_r;
1106 cs_real_t theta_met_min= cs_math_big_r;
1107
1108 cs_real_t *z_ground = NULL;
1109 if (aopt->compute_z_ground == true) {
1110 cs_atmo_z_ground_compute();
1111 z_ground = cs_field_by_name_try("z_ground")->val;
1112 }
1113
1114 BFT_MALLOC(dlmo_var, m->n_cells_with_ghosts, cs_real_t);
1115
1116 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells_with_ghosts; cell_id++) {
1117 dlmo_var[cell_id] = 0.0;
1118 }
1119 /* For DRSM models store Rxz/k */
1120 cs_field_t *f_axz = cs_field_by_name_try("meteo_shear_anisotropy");
1121
1122 if (dlmo > 0)
1123 ri_max = 0.75; // Value chosen to limit buoyancy vs shear production
1124
1125 /* Profiles */
1126 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1127
1128 cs_real_t z_grd = 0.;
1129 if (z_ground != NULL)
1130 z_grd = z_ground[cell_id];
1131
1132 /* Local elevation */
1133 cs_real_t z = cell_cen[cell_id][2] - z_grd;
1134
1135 /* Velocity profile */
1136 cs_real_t u_norm = ustar0 / kappa * cs_mo_psim(z+z0, z0, dlmo);
1137
1138 cpro_met_vel[cell_id][0] = - sin(angle * cs_math_pi/180.) * u_norm;
1139 cpro_met_vel[cell_id][1] = - cos(angle * cs_math_pi/180.) * u_norm;
1140
1141 /* Potential temperature profile
1142 * Note: same roughness as dynamics */
1143 cpro_met_potemp[cell_id] = theta0 + tstar / kappa * cs_mo_psih(z+z0, z0, dlmo);
1144
1145 /* Richardson flux number profile */
1146 // Note : ri_f = z/(Pr_t L) * phih/phim^2 = z/Lmo / phim
1147 cs_real_t ri_f = (z+z0) * dlmo / cs_mo_phim(z+z0, dlmo);
1148
1149 /* TKE profile */
1150 cpro_met_k[cell_id] = cs_math_pow2(ustar0) / sqrt(cmu)
1151 * sqrt(1. - CS_MIN(ri_f, 1.));
1152 if (f_axz != NULL)
1153 f_axz->val[cell_id] = -sqrt(cmu / (1. - CS_MIN(ri_f, ri_max)));
1154
1155 /* epsilon profile */
1156 cpro_met_eps[cell_id] = cs_math_pow3(ustar0) / (kappa * (z+z0))
1157 * cs_mo_phim(z+z0, dlmo)*(1.-CS_MIN(ri_f, 1.));
1158
1159 /* Very stable cases */
1160 if (ri_f > ri_max) {
1161 if (z < z_lim) {
1162 //Ri_f is an increasing monotonic function, so the lowest value of
1163 //z for which Ri_f>Ri_max is needed
1164 z_lim = z;
1165 u_met_min=u_norm;
1166 theta_met_min=cpro_met_potemp[cell_id];
1167 }
1168 }
1169 }
1170
1171 cs_parall_min(1,CS_REAL_TYPE, &z_lim);
1172 cs_parall_min(1,CS_REAL_TYPE, &u_met_min);
1173 cs_parall_min(1,CS_REAL_TYPE, &theta_met_min);
1174
1175 /* Very stable cases, corresponding to mode 0 in the Python prepro */
1176 if (z_lim < 0.5*cs_math_big_r) { // Clipping only if there are cells to be clipped
1177 bft_printf("Switching to very stable clipping for meteo profile.\n");
1178 bft_printf("All altitudes above %f have been modified by clipping.\n",z_lim);
1179 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1180
1181 cs_real_t z_grd = 0.;
1182 if (z_ground != NULL)
1183 z_grd = z_ground[cell_id];
1184
1185 cs_real_t z = cell_cen[cell_id][2] - z_grd;
1186 if (z >= z_lim) {
1187 /* mode = 0 is ustar=cst */
1188 dlmo_var[cell_id] = dlmo * (z_lim + z0) / (z + z0);
1189
1190 /* Velocity profile */
1191 cs_real_t u_norm = u_met_min + ustar0 / kappa
1192 * cs_mo_phim(z_lim + z0, dlmo)
1193 * log((z+z0) / (z_lim+z0));
1194
1195 cpro_met_vel[cell_id][0] = - sin(angle * cs_math_pi/180.) * u_norm;
1196 cpro_met_vel[cell_id][1] = - cos(angle * cs_math_pi/180.) * u_norm;
1197
1198 /* Potential temperature profile
1199 * Note: same roughness as dynamics */
1200 cpro_met_potemp[cell_id] = theta_met_min
1201 + tstar * (z_lim+z0) / kappa
1202 * cs_mo_phih(z_lim+z0, dlmo)
1203 * (-1./(z+z0) + 1./(z_lim+z0)) ;
1204 /* TKE profile
1205 ri_max is necessarily lower than 1, but CS_MIN might be useful if
1206 that changes in the future */
1207 cpro_met_k[cell_id] = cs_math_pow2(ustar0) / sqrt(cmu)
1208 * sqrt(1. - CS_MIN(ri_max, 1.));
1209
1210 if (f_axz != NULL)
1211 f_axz->val[cell_id] = -sqrt(cmu / (1. - CS_MIN(ri_max, 1.)));
1212
1213 /* epsilon profile */
1214 cpro_met_eps[cell_id] = cs_math_pow3(ustar0) / kappa * dlmo_var[cell_id]
1215 * (1- CS_MIN(ri_max, 1.)) / CS_MIN(ri_max, 1.);
1216 }
1217 }
1218 }
1219
1220 cs_atmo_hydrostatic_profiles_compute();
1221 BFT_FREE(dlmo_var);
1222 }
1223
1224 /*----------------------------------------------------------------------------*/
1225 /*!
1226 * \brief Compute the ground elevation.
1227 *
1228 * This function solves the following transport equation on \f$ \varia \f$:
1229 * \f[
1230 * \dfrac{\partial \varia}{\partial t} + \divs \left( \varia \vect{g} \right)
1231 * - \divs \left( \vect{V} \right) \varia = 0
1232 * \f]
1233 * where \f$ \vect{g} \f$ is the gravity field
1234 *
1235 * The boundary conditions on \f$ \varia \f$ read:
1236 * \f[
1237 * \varia = z \textrm{ on walls}
1238 * \f]
1239 * \f[
1240 * \dfrac{\partial \varia}{\partial n} = 0 \textrm{ elsewhere}
1241 * \f]
1242 *
1243 * Remarks:
1244 * - a steady state is looked for.
1245 */
1246 /*----------------------------------------------------------------------------*/
1247
1248 void
cs_atmo_z_ground_compute(void)1249 cs_atmo_z_ground_compute(void)
1250 {
1251 /* Initialization
1252 *===============*/
1253
1254 if (!_atmo_option.compute_z_ground)
1255 return;
1256
1257 const cs_domain_t *domain = cs_glob_domain;
1258 const cs_mesh_t *m = domain->mesh;
1259 const cs_mesh_quantities_t *mq = domain->mesh_quantities;
1260
1261 const cs_real_3_t *restrict i_face_normal =
1262 (const cs_real_3_t *restrict)mq->i_face_normal;
1263 const cs_real_3_t *restrict b_face_normal =
1264 (const cs_real_3_t *restrict)mq->b_face_normal;
1265 const cs_real_3_t *restrict b_face_cog
1266 = (const cs_real_3_t *restrict)mq->b_face_cog;
1267
1268 const int *bc_type = cs_glob_bc_type;
1269
1270 /* Pointer to z_ground field */
1271 cs_field_t *f = cs_field_by_name_try("z_ground");
1272
1273 cs_real_t *restrict i_massflux
1274 = cs_field_by_id
1275 (cs_field_get_key_int(f, cs_field_key_id("inner_mass_flux_id")))->val;
1276 cs_real_t *restrict b_massflux
1277 = cs_field_by_id
1278 (cs_field_get_key_int(f, cs_field_key_id("boundary_mass_flux_id")))->val;
1279
1280 cs_var_cal_opt_t vcopt;
1281 cs_field_get_key_struct(f, cs_field_key_id("var_cal_opt"), &vcopt);
1282
1283 cs_real_t normal[3];
1284 /* Normal direction is given by the gravity */
1285 cs_math_3_normalise((const cs_real_t *)(cs_glob_physical_constants->gravity),
1286 normal);
1287
1288 for (int i = 0; i < 3; i++)
1289 normal[i] *= -1.;
1290
1291 /* Compute the mass flux due to V = - g / ||g||
1292 *=============================================*/
1293
1294 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++)
1295 i_massflux[face_id] = cs_math_3_dot_product(normal, i_face_normal[face_id]);
1296
1297 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++)
1298 b_massflux[face_id] = cs_math_3_dot_product(normal, b_face_normal[face_id]);
1299
1300 /* Boundary conditions
1301 *====================*/
1302
1303 cs_real_t norm = 0.;
1304 cs_real_t ground_surf = 0.;
1305
1306 /* Dirichlet at walls, homogeneous Neumann elsewhere */
1307
1308 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1309 /* Dirichlet BCs */
1310 if ((bc_type[face_id] == CS_SMOOTHWALL || bc_type[face_id] == CS_ROUGHWALL)
1311 && b_massflux[face_id] <= 0.) {
1312
1313 vcopt.ndircl = 1;
1314 cs_real_t hint = 1. / mq->b_dist[face_id];
1315 cs_real_t pimp = cs_math_3_dot_product(b_face_cog[face_id], normal);
1316
1317 cs_boundary_conditions_set_dirichlet_scalar(&(f->bc_coeffs->a[face_id]),
1318 &(f->bc_coeffs->af[face_id]),
1319 &(f->bc_coeffs->b[face_id]),
1320 &(f->bc_coeffs->bf[face_id]),
1321 pimp,
1322 hint,
1323 cs_math_infinite_r);
1324 norm += cs_math_pow2(f->bc_coeffs->a[face_id]) * mq->b_face_surf[face_id];
1325 ground_surf += mq->b_face_surf[face_id];
1326 }
1327 /* Neumann Boundary Conditions */
1328 else {
1329
1330 cs_real_t hint = 1. / mq->b_dist[face_id];
1331 cs_real_t qimp = 0.;
1332
1333 cs_boundary_conditions_set_neumann_scalar(&(f->bc_coeffs->a[face_id]),
1334 &(f->bc_coeffs->af[face_id]),
1335 &(f->bc_coeffs->b[face_id]),
1336 &(f->bc_coeffs->bf[face_id]),
1337 qimp,
1338 hint);
1339
1340 }
1341 }
1342
1343 cs_parall_max(1, CS_INT_TYPE, &(vcopt.ndircl));
1344
1345 /* Matrix
1346 *=======*/
1347
1348 cs_real_t *rovsdt, *dpvar;
1349 BFT_MALLOC(rovsdt, m->n_cells_with_ghosts, cs_real_t);
1350 BFT_MALLOC(dpvar, m->n_cells_with_ghosts, cs_real_t);
1351
1352 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells_with_ghosts; cell_id++) {
1353 rovsdt[cell_id] = 0.;
1354 dpvar[cell_id] = 0.;
1355 }
1356 /* Right hand side
1357 *================*/
1358
1359 cs_real_t *rhs;
1360 BFT_MALLOC(rhs, m->n_cells_with_ghosts, cs_real_t);
1361
1362 for (cs_lnum_t cell_id = 0; cell_id< m->n_cells_with_ghosts; cell_id++)
1363 rhs[cell_id] = 0.;
1364
1365 /* Norm
1366 *======*/
1367
1368 cs_parall_sum(1, CS_REAL_TYPE, &norm);
1369 cs_parall_sum(1, CS_REAL_TYPE, &ground_surf);
1370
1371 if (ground_surf > 0.)
1372 norm = sqrt(norm / ground_surf) * mq->tot_vol;
1373 else {
1374 bft_printf("No ground BC or no gravity: no computation of ground elevation.\n");
1375 return;
1376 }
1377
1378 /* Solving
1379 *=========*/
1380
1381 /* In case of a theta-scheme, set theta = 1;
1382 no relaxation in steady case either */
1383 cs_real_t inf_norm = 1.;
1384
1385 /* Overall loop in order to ensure convergence */
1386 for (int sweep = 0; sweep < vcopt.nswrsm && inf_norm > vcopt.epsrsm; sweep++) {
1387
1388 cs_equation_iterative_solve_scalar(0, /* idtvar: no steady state algo */
1389 -1, /* no over loops */
1390 f->id,
1391 f->name,
1392 0, /* iescap */
1393 0, /* imucpp */
1394 norm,
1395 &vcopt,
1396 f->val_pre,
1397 f->val,
1398 f->bc_coeffs->a,
1399 f->bc_coeffs->b,
1400 f->bc_coeffs->af,
1401 f->bc_coeffs->bf,
1402 i_massflux,
1403 b_massflux,
1404 i_massflux, /* viscosity, not used */
1405 b_massflux, /* viscosity, not used */
1406 i_massflux, /* viscosity, not used */
1407 b_massflux, /* viscosity, not used */
1408 NULL,
1409 NULL,
1410 NULL,
1411 0, /* icvflb (upwind) */
1412 NULL,
1413 rovsdt,
1414 rhs,
1415 f->val,
1416 dpvar,
1417 NULL,
1418 NULL);
1419
1420 /* Compute the L_infinity norm */
1421 inf_norm = 0.;
1422 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1423 //FIXME make this dimensionless
1424 inf_norm = fmax(inf_norm, fabs(f->val[cell_id]-f->val_pre[cell_id]));
1425
1426 /* Current to previous */
1427 f->val_pre[cell_id] = f->val[cell_id];
1428 }
1429
1430 cs_parall_max(1, CS_REAL_TYPE, &inf_norm);
1431 }
1432
1433 /* Free memory */
1434 BFT_FREE(dpvar);
1435 BFT_FREE(rhs);
1436 BFT_FREE(rovsdt);
1437 }
1438
1439 /*----------------------------------------------------------------------------*/
1440 /*!
1441 * \brief Compute hydrostatic profiles of density and pressure.
1442 *
1443 * This function solves the following transport equation on \f$ \varia \f$:
1444 * \f[
1445 * \divs \left( \grad \varia \right)
1446 * = \divs \left( \dfrac{\vect{g}}{c_p \theta} \right)
1447 * \f]
1448 * where \f$ \vect{g} \f$ is the gravity field and \f$ \theta \f$
1449 * is the potential temperature.
1450 *
1451 * The boundary conditions on \f$ \varia \f$ read:
1452 * \f[
1453 * \varia = \left(\dfrac{P_{sea}}{p_s}\right)^{R/C_p} \textrm{on the ground}
1454 * \f]
1455 * and Neumann elsewhere.
1456 */
1457 /*----------------------------------------------------------------------------*/
1458
1459 void
cs_atmo_hydrostatic_profiles_compute(void)1460 cs_atmo_hydrostatic_profiles_compute(void)
1461 {
1462 /* Initialization
1463 *===============*/
1464
1465 cs_mesh_t *m = cs_glob_mesh;
1466 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1467
1468 const cs_real_3_t *restrict b_face_cog
1469 = (const cs_real_3_t *restrict)mq->b_face_cog;
1470 const cs_real_3_t *restrict cell_cen
1471 = (const cs_real_3_t *restrict)mq->cell_cen;
1472
1473 cs_physical_constants_t *phys_cst = cs_get_glob_physical_constants();
1474 cs_field_t *f = cs_field_by_name("meteo_pressure");
1475 cs_var_cal_opt_t vcopt;
1476 cs_field_get_key_struct(f, cs_field_key_id("var_cal_opt"), &vcopt);
1477 cs_field_t *potemp = cs_field_by_name("meteo_pot_temperature");
1478 cs_field_t *density = cs_field_by_name("meteo_density");
1479 cs_field_t *temp = cs_field_by_name("meteo_temperature");
1480
1481 cs_atmo_option_t *aopt = &_atmo_option;
1482 cs_real_t g = cs_math_3_norm(phys_cst->gravity);
1483
1484 const cs_fluid_properties_t *phys_pro = cs_get_glob_fluid_properties();
1485 /* potential temp at ref */
1486 cs_real_t pref = cs_glob_atmo_constants->ps;
1487 cs_real_t rair = phys_pro->r_pg_cnst;
1488 cs_real_t cp0 = phys_pro->cp0;
1489 cs_real_t rscp = rair/cp0; /* Around 2/7 */
1490 cs_real_t theta0 = aopt->meteo_t0 * pow(pref/ aopt->meteo_psea, rscp);
1491
1492 const cs_velocity_pressure_model_t *vp_model
1493 = cs_glob_velocity_pressure_model;
1494 const int idilat = vp_model->idilat;
1495
1496 /* Get the lowest altitude (should also be minimum of z_ground)
1497 *=============================================================*/
1498
1499 cs_real_t z_min = cs_math_big_r;
1500
1501 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id ++) {
1502 if (b_face_cog[face_id][2] < z_min)
1503 z_min = b_face_cog[face_id][2];
1504 }
1505
1506 cs_parall_min(1, CS_REAL_TYPE, &z_min);
1507
1508 /* p_ground is pressure at the lowest level */
1509
1510 cs_real_t p_ground = aopt->meteo_psea;
1511
1512 /* Initialize temperature, pressure and density from neutral conditions
1513 *=====================================================================*/
1514 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1515 cs_real_t z = cell_cen[cell_id][2] - z_min;
1516 cs_real_t zt = fmin(z, 11000.);
1517 cs_real_t factor = fmax(1. - g * zt / (cp0 * aopt->meteo_t0), 0.);
1518 temp->val[cell_id] = aopt->meteo_t0 * factor ;
1519 f->val[cell_id] = p_ground * pow(factor, rscp)
1520 * exp(- g/(rair*temp->val[cell_id]) * (z - zt));
1521 density->val[cell_id] = f->val[cell_id] / (rair * temp->val[cell_id]);
1522 }
1523
1524 /* Boussinesq hypothesis */
1525 if (idilat==0) {
1526 bft_printf(
1527 "Meteo profiles are computed according to Boussinesq approximation.\n"
1528 "Using adiabatic profiles for temperature and pressure."
1529 "Density is computed accordingly.\n");
1530 }
1531
1532 cs_real_t *i_massflux = NULL;
1533 cs_real_t *b_massflux = NULL;
1534 BFT_MALLOC(i_massflux, m->n_i_faces, cs_real_t);
1535 BFT_MALLOC(b_massflux, m->n_b_faces, cs_real_t);
1536
1537 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++)
1538 i_massflux[face_id] = 0.;
1539
1540 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++)
1541 b_massflux[face_id] = 0.;
1542
1543 cs_real_t *i_viscm = NULL;
1544 BFT_MALLOC(i_viscm, m->n_i_faces, cs_real_t);
1545
1546 cs_real_t *b_viscm = NULL;
1547 BFT_MALLOC(b_viscm, m->n_b_faces, cs_real_t);
1548
1549 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++)
1550 i_viscm[face_id] = 0.;
1551
1552 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++)
1553 b_viscm[face_id] = 0.;
1554
1555 cs_real_3_t *f_ext, *dfext;
1556 BFT_MALLOC(f_ext, m->n_cells_with_ghosts, cs_real_3_t);
1557 BFT_MALLOC(dfext, m->n_cells_with_ghosts, cs_real_3_t);
1558
1559 /* dfext is actually a dummy used to copy calhyd */
1560 /* f_ext is initialized with an initial density */
1561
1562 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1563 f_ext[cell_id][0] = density->val[cell_id] * phys_cst->gravity[0];
1564 dfext[cell_id][0] = 0.;
1565 f_ext[cell_id][1] = density->val[cell_id] * phys_cst->gravity[1];
1566 dfext[cell_id][1] = 0.;
1567 f_ext[cell_id][2] = density->val[cell_id] * phys_cst->gravity[2];
1568 dfext[cell_id][2] = 0;
1569 }
1570
1571 /* Solving
1572 *=========*/
1573
1574 cs_real_t *dam = NULL;
1575 BFT_MALLOC(dam, m->n_cells_with_ghosts, cs_real_t);
1576
1577 cs_real_t *xam = NULL;
1578 BFT_MALLOC(xam, m->n_i_faces, cs_real_t);
1579
1580 cs_real_t *rhs = NULL;
1581 BFT_MALLOC(rhs, m->n_cells_with_ghosts, cs_real_t);
1582
1583 cs_real_t *dpvar = NULL;
1584 BFT_MALLOC(dpvar, m->n_cells_with_ghosts, cs_real_t);
1585
1586 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1587 dpvar[cell_id] = 0.;
1588 dam[cell_id] = 0.;
1589 rhs[cell_id] = 0.;
1590 }
1591
1592 cs_real_t inf_norm = 1.;
1593
1594 /* Loop to compute pressure profile */
1595 for (int sweep = 0; sweep < vcopt.nswrsm && inf_norm > vcopt.epsrsm; sweep++) {
1596 //FIXME 100 or nswrsm
1597
1598 /* Update previous values of pressure for the convergence test */
1599 cs_field_current_to_previous(f);
1600
1601 _hydrostatic_pressure_compute(f_ext,
1602 dfext,
1603 f->val,
1604 f->name,
1605 z_min,
1606 p_ground,
1607 i_massflux,
1608 b_massflux,
1609 i_viscm,
1610 b_viscm,
1611 dam,
1612 xam,
1613 dpvar,
1614 rhs);
1615
1616 /* L infinity residual computation and forcing update */
1617 inf_norm = 0.;
1618 for (cs_lnum_t cell_id = 0; cell_id < m->n_cells; cell_id++) {
1619 inf_norm = fmax(fabs(f->val[cell_id] - f->val_pre[cell_id])/pref, inf_norm);
1620
1621 /* Boussinesq hypothesis: do not update adiabatic temperature profile */
1622 if (idilat > 0) {
1623 temp->val[cell_id] = potemp->val[cell_id]
1624 * pow((f->val[cell_id]/pref), rscp);
1625 }
1626
1627 /* f_ext = rho^k * g */
1628 cs_real_t rho_k = f->val[cell_id] / (rair * temp->val[cell_id]);
1629 density->val[cell_id] = rho_k;
1630
1631 f_ext[cell_id][0] = rho_k * phys_cst->gravity[0];
1632 f_ext[cell_id][1] = rho_k * phys_cst->gravity[1];
1633 f_ext[cell_id][2] = rho_k * phys_cst->gravity[2];
1634 }
1635 cs_parall_max(1, CS_REAL_TYPE ,&inf_norm);
1636
1637 if (cs_log_default_is_active())
1638 bft_printf
1639 (_("Meteo profiles: iterative process to compute hydrostatic pressure\n"
1640 " sweep %d, L infinity norm (delta p) / ps =%e\n"), sweep, inf_norm);
1641 }
1642
1643 /* Free memory */
1644 BFT_FREE(dpvar);
1645 BFT_FREE(rhs);
1646 BFT_FREE(i_viscm);
1647 BFT_FREE(b_viscm);
1648 BFT_FREE(xam);
1649 BFT_FREE(dam);
1650 BFT_FREE(f_ext);
1651 BFT_FREE(dfext);
1652 BFT_FREE(b_massflux);
1653 BFT_FREE(i_massflux);
1654 }
1655
1656 /*----------------------------------------------------------------------------*/
1657 /*!
1658 * \brief This function set the file name of the meteo file.
1659 *
1660 * \param[in] file_name name of the file.
1661 */
1662 /*----------------------------------------------------------------------------*/
1663
1664 void
cs_atmo_set_meteo_file_name(const char * file_name)1665 cs_atmo_set_meteo_file_name(const char *file_name)
1666 {
1667 if (file_name == NULL) {
1668 return;
1669 }
1670
1671 if (_atmo_option.meteo_file_name == NULL) {
1672 BFT_MALLOC(_atmo_option.meteo_file_name,
1673 strlen(file_name) + 1,
1674 char);
1675 }
1676 else {
1677 BFT_REALLOC(_atmo_option.meteo_file_name,
1678 strlen(file_name) + 1,
1679 char);
1680 }
1681
1682 sprintf(_atmo_option.meteo_file_name, "%s", file_name);
1683 }
1684
1685 /*----------------------------------------------------------------------------*/
1686 /*!
1687 * \brief This function set the file name of the chemistry file.
1688 *
1689 * \param[in] file_name name of the file.
1690 */
1691 /*----------------------------------------------------------------------------*/
1692
1693 void
cs_atmo_set_chem_conc_file_name(const char * file_name)1694 cs_atmo_set_chem_conc_file_name(const char *file_name)
1695 {
1696 if (file_name == NULL) {
1697 return;
1698 }
1699
1700 if (_atmo_chem.chem_conc_file_name == NULL) {
1701 BFT_MALLOC(_atmo_chem.chem_conc_file_name,
1702 strlen(file_name) + 1,
1703 char);
1704 }
1705 else {
1706 BFT_REALLOC(_atmo_chem.chem_conc_file_name,
1707 strlen(file_name) + 1,
1708 char);
1709 }
1710
1711 sprintf(_atmo_chem.chem_conc_file_name, "%s", file_name);
1712 }
1713
1714 /*----------------------------------------------------------------------------*/
1715 /*!
1716 * \brief This function set the file name of the aerosol file.
1717 *
1718 * \param[in] file_name name of the file.
1719 */
1720 /*----------------------------------------------------------------------------*/
1721
1722 void
cs_atmo_set_aero_conc_file_name(const char * file_name)1723 cs_atmo_set_aero_conc_file_name(const char *file_name)
1724 {
1725 if (file_name == NULL) {
1726 return;
1727 }
1728 if (cs_glob_atmo_chemistry->aerosol_model == CS_ATMO_AEROSOL_OFF) {
1729 return;
1730 }
1731
1732 if (_atmo_chem.aero_conc_file_name == NULL) {
1733 BFT_MALLOC(_atmo_chem.aero_conc_file_name,
1734 strlen(file_name) + 1,
1735 char);
1736 }
1737 else {
1738 BFT_REALLOC(_atmo_chem.aero_conc_file_name,
1739 strlen(file_name) + 1,
1740 char);
1741 }
1742
1743 sprintf(_atmo_chem.aero_conc_file_name, "%s", file_name);
1744 }
1745
1746 /*----------------------------------------------------------------------------*/
1747 /*!
1748 * \brief This function set the file name of the SPACK file.
1749 *
1750 * \param[in] file_name name of the file.
1751 */
1752 /*----------------------------------------------------------------------------*/
1753
1754 void
cs_atmo_chemistry_set_spack_file_name(const char * file_name)1755 cs_atmo_chemistry_set_spack_file_name(const char *file_name)
1756 {
1757 if (file_name == NULL) {
1758 _atmo_chem.model = 0;
1759 return;
1760 }
1761
1762 _atmo_chem.model = 4;
1763
1764 BFT_MALLOC(_atmo_chem.spack_file_name,
1765 strlen(file_name) + 1,
1766 char);
1767
1768 sprintf(_atmo_chem.spack_file_name, "%s", file_name);
1769 }
1770
1771 /*----------------------------------------------------------------------------*/
1772 /*!
1773 * \brief This function sets the file name to initialize the aerosol library.
1774 *
1775 * \param[in] file_name name of the file.
1776 */
1777 /*----------------------------------------------------------------------------*/
1778
1779 void
cs_atmo_chemistry_set_aerosol_file_name(const char * file_name)1780 cs_atmo_chemistry_set_aerosol_file_name(const char *file_name)
1781 {
1782 if (file_name == NULL) {
1783 _atmo_chem.aerosol_model = CS_ATMO_AEROSOL_OFF;
1784 return;
1785 }
1786
1787 _atmo_chem.aerosol_model = CS_ATMO_AEROSOL_SSH;
1788
1789 BFT_MALLOC(_atmo_chem.aero_file_name,
1790 strlen(file_name) + 1,
1791 char);
1792
1793 sprintf(_atmo_chem.aero_file_name, "%s", file_name);
1794 }
1795
1796 /*----------------------------------------------------------------------------*/
1797 /*!
1798 * \brief This function declares additional transported variables for
1799 * atmospheric module for the chemistry defined from SPACK.
1800 */
1801 /*----------------------------------------------------------------------------*/
1802
1803 void
cs_atmo_declare_chem_from_spack(void)1804 cs_atmo_declare_chem_from_spack(void)
1805 {
1806 assert(_atmo_chem.model == 4);
1807
1808 if (_atmo_chem.spack_file_name == NULL)
1809 bft_error(__FILE__,__LINE__, 0,
1810 _("Atmo chemistry from SPACK file: requires a SPACK file."));
1811
1812 char line[512] = "";
1813
1814 /* Open file */
1815 bft_printf("SPACK file for atmo chemistry:\n %s \n",
1816 _atmo_chem.spack_file_name);
1817
1818 FILE *file = fopen(_atmo_chem.spack_file_name, "rt");
1819 if (file == NULL)
1820 bft_error(__FILE__,__LINE__, 0,
1821 _("Atmo chemistry from SPACK file: Could not open file."));
1822
1823 int line_num = 0;
1824
1825 /* Read "[species]" */
1826 while (true) {
1827 line_num++;
1828 if (fscanf(file, "%s\n", line) != 1)
1829 bft_error(__FILE__,__LINE__, 0,
1830 _("Atmo chemistry from SPACK file: Could not skip header."));
1831
1832 if (strncmp("[species]", line, 512) == 0)
1833 break;
1834 }
1835
1836 /* Read SPACK: first loop count the number of species */
1837 for (int i = 1; true; i++ ) {
1838 /* Read species */
1839 line_num++;
1840 if (fscanf(file, "%s\n", line) != 1)
1841 bft_error(__FILE__,__LINE__, 0,
1842 _("Atmo chemistry from SPACK file: Could not read line %d."),
1843 line_num);
1844
1845 /* When reach [molecular_waight]: break */
1846 if (strncmp("[molecular_weight]", line, 512) == 0)
1847 break;
1848 else {
1849 _atmo_chem.n_species = i;
1850 }
1851 }
1852
1853 /* Now allocate arrays */
1854 BFT_MALLOC(_atmo_chem.species_to_field_id, _atmo_chem.n_species, int);
1855 BFT_MALLOC(_atmo_chem.species_to_scalar_id, _atmo_chem.n_species, int);
1856 BFT_MALLOC(_atmo_chem.molar_mass, _atmo_chem.n_species, cs_real_t);
1857 BFT_MALLOC(_atmo_chem.chempoint, _atmo_chem.n_species, int);
1858
1859 /* Read SPACK: second loop Create variables and read molar mass */
1860 for (int i = 0; i < _atmo_chem.n_species; i++ ) {
1861 char name[512] = "";
1862 char label[512] = "";
1863
1864 /* Read species name and molar mass */
1865 line_num++;
1866 if (fscanf(file, "%s %lf\n", line, &(_atmo_chem.molar_mass[i])) != 2)
1867 bft_error(__FILE__,__LINE__, 0,
1868 _("Atmospheric chemistry from SPACK file:\n"
1869 " could not read species name and molar mass, line %d."),
1870 line_num);
1871
1872 /* The order is already ok */
1873 _atmo_chem.chempoint[i] = i+1; //FIXME ?
1874
1875 /* Build name of the field:
1876 * species_name in lower case */
1877 strcpy(name, "species_");
1878 _strtolower(label, line);
1879 strcat(name, label);
1880
1881 /* Field of dimension 1 */
1882 /* Give the original name as label */
1883 _atmo_chem.species_to_field_id[i]
1884 = cs_variable_field_create(name, line, CS_MESH_LOCATION_CELLS, 1);
1885
1886 /* Scalar field, store in isca_chem/species_to_scalar_id (FORTRAN/C) array */
1887 _atmo_chem.species_to_scalar_id[i]
1888 = cs_add_model_field_indexes(_atmo_chem.species_to_field_id[i]);
1889
1890 }
1891 }
1892
1893 /*----------------------------------------------------------------------------*/
1894 /*!
1895 * \brief 1D Radiative scheme - Solar data + zenithal angle)
1896 *
1897 * Compute:
1898 * - zenithal angle
1899 * - solar contant (with correction for earth - solar length)
1900 * - albedo if above the sea
1901 * (Use analytical formulae of Paltrige and Platt
1902 * dev.in atm. science no 5)
1903 * \param[in] latitude latitude
1904 * \param[in] longitude longitude
1905 * \param[in] squant start day in the year
1906 * \param[in] utc Universal time (hour)
1907 * \param[in] sea_id sea index
1908 * \param[out] albedo albedo
1909 * \param[out] muzero cosin of zenithal angle
1910 * \param[out] omega solar azimut angle
1911 * \param[out] fo solar constant
1912 */
1913 /*----------------------------------------------------------------------------*/
1914
1915 void
cs_atmo_compute_solar_angles(cs_real_t latitude,cs_real_t longitude,cs_real_t squant,cs_real_t utc,int sea_id,cs_real_t * albedo,cs_real_t * muzero,cs_real_t * omega,cs_real_t * fo)1916 cs_atmo_compute_solar_angles(cs_real_t latitude,
1917 cs_real_t longitude,
1918 cs_real_t squant,
1919 cs_real_t utc,
1920 int sea_id,
1921 cs_real_t *albedo,
1922 cs_real_t *muzero,
1923 cs_real_t *omega,
1924 cs_real_t *fo)
1925 {
1926 /* 1 - initialisations */
1927 *fo = 1370.;
1928
1929 /* conversions sexagesimal-decimal */
1930
1931 cs_real_t flat = latitude *cs_math_pi/180.;
1932 cs_real_t flong = longitude * 4. / 60.;
1933
1934 cs_real_t t00 = 2. * cs_math_pi * squant/365.;
1935
1936 /* 2 - compute declinaison (maximum error < 3 mn) */
1937
1938 cs_real_t decl
1939 = 0.006918 - 0.399912*cos(t00) + 0.070257*sin(t00)
1940 - 0.006758*cos(2.*t00) + 0.000907*sin(2.*t00) - 0.002697*cos(3.*t00)
1941 + 0.001480*sin(3.*t00);
1942
1943 /* 3 - compute local solar hour
1944 * equation du temps erreur maxi 35 secondes
1945 */
1946
1947 cs_real_t eqt
1948 = ( 0.000075 + 0.001868*cos(t00) - 0.032077*sin(t00)
1949 - 0.014615*cos(2.*t00) - 0.040849*sin(2.*t00)) * 12./cs_math_pi;
1950
1951 cs_real_t local_time = utc + flong + eqt;
1952
1953 /* Transformation local_time-radians */
1954
1955 /* On retire cs_math_pi et on prend le modulo 2pi du resultat */
1956 cs_real_t hr = (local_time - 12.)*cs_math_pi/12.;
1957 if (local_time < 12.)
1958 hr = (local_time + 12.)*cs_math_pi/12.;
1959
1960 /* 4 - compute of cosinus of the zenitghal angle */
1961
1962 *muzero = sin(decl)*sin(flat) + cos(decl)*cos(flat)*cos(hr);
1963
1964 cs_real_t za = acos(*muzero);
1965
1966 /* 5 - compute solar azimut */
1967 *omega = 0.;
1968 if (CS_ABS(sin(za)) > cs_math_epzero) {
1969 /* Cosinus of the zimut angle */
1970 cs_real_t co = (sin(decl)*cos(flat)-cos(decl)*sin(flat)*cos(hr))/sin(za);
1971 *omega = acos(co);
1972 if (local_time > 12.)
1973 *omega = 2. * cs_math_pi - acos(co);
1974 }
1975 *omega -= cs_glob_atmo_option->domain_orientation * cs_math_pi / 180.;
1976
1977 /* 5 - compute albedo at sea which depends on the zenithal angle */
1978
1979 if (sea_id == 1) {
1980 cs_real_t ho = acos(*muzero);
1981 ho = 180.*(cs_math_pi/2. - ho)/cs_math_pi;
1982 if (ho < 8.5)
1983 ho = 8.5;
1984 if (ho > 60.)
1985 ho = 60.;
1986 *albedo = 3./ho;
1987 }
1988
1989 /* 6 - Compute solar constant
1990 distance correction earth-sun
1991 corfo=(r0/r)**2
1992 precision better than e-04 */
1993
1994 cs_real_t corfo
1995 = 1.000110 + 0.034221*cos(t00) + 0.001280*sin(t00)
1996 + 0.000719*cos(2.*t00) + 0.000077*sin(2.*t00);
1997 *fo *= corfo;
1998 }
1999
2000 /*----------------------------------------------------------------------------*/
2001 /*!
2002 * \brief Print the atmospheric module options to setup.log.
2003 */
2004 /*----------------------------------------------------------------------------*/
2005
2006 void
cs_atmo_log_setup(void)2007 cs_atmo_log_setup(void)
2008 {
2009 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] == CS_ATMO_OFF)
2010 return;
2011
2012 cs_log_printf(CS_LOG_SETUP,
2013 _("\n"
2014 "Atmospheric module options\n"
2015 "--------------------------\n\n"));
2016
2017 switch(cs_glob_physical_model_flag[CS_ATMOSPHERIC]) {
2018 case CS_ATMO_CONSTANT_DENSITY:
2019 /* Constant density */
2020 cs_log_printf(CS_LOG_SETUP,
2021 _(" Constant density\n\n"));
2022 break;
2023 case CS_ATMO_DRY:
2024 /* Dry */
2025 cs_log_printf(CS_LOG_SETUP,
2026 _(" Dry atmosphere\n\n"));
2027 break;
2028 case CS_ATMO_HUMID:
2029 /* Humid */
2030 cs_log_printf(CS_LOG_SETUP,
2031 _(" Humid atmosphere\n\n"));
2032 break;
2033 default:
2034 break;
2035 }
2036
2037 if (cs_glob_atmo_option->compute_z_ground > 0)
2038 cs_log_printf(CS_LOG_SETUP,
2039 _(" Compute ground elevation\n\n"));
2040
2041 if (cs_glob_atmo_option->open_bcs_treatment > 0)
2042 cs_log_printf(CS_LOG_SETUP,
2043 _(" Impose open BCs with momentum source terms\n"));
2044
2045 if (cs_glob_atmo_option->open_bcs_treatment == 2)
2046 cs_log_printf(CS_LOG_SETUP,
2047 _(" and impose profiles at ingoing faces\n\n"));
2048
2049 /* CUT */
2050 cs_log_printf
2051 (CS_LOG_SETUP,
2052 _("\n"
2053 " Starting Coordinated Universal Time (CUT):\n"
2054 " Year: %4d\n"
2055 " Quant: %4d\n"
2056 " Hour: %4d\n"
2057 " Min: %4d\n"
2058 " Sec: %4f\n\n"),
2059 cs_glob_atmo_option->syear,
2060 cs_glob_atmo_option->squant,
2061 cs_glob_atmo_option->shour,
2062 cs_glob_atmo_option->smin,
2063 cs_glob_atmo_option->ssec);
2064
2065 /* Centre of the domain latitude */
2066 cs_log_printf
2067 (CS_LOG_SETUP,
2068 _(" Domain center:\n"
2069 " latitude: %6f\n"
2070 " longitude: %6f\n"
2071 " x center (in Lambert-93) : %6f\n"
2072 " y center (in Lmabert-93) : %6f\n\n"),
2073 cs_glob_atmo_option->latitude,
2074 cs_glob_atmo_option->longitude,
2075 cs_glob_atmo_option->x_l93,
2076 cs_glob_atmo_option->y_l93);
2077
2078 if (cs_glob_atmo_option->meteo_profile == 1) {
2079 cs_log_printf
2080 (CS_LOG_SETUP,
2081 _(" Large scale Meteo file: %s\n\n"),
2082 cs_glob_atmo_option->meteo_file_name);
2083 }
2084
2085
2086 if (cs_glob_atmo_option->meteo_profile == 2) {
2087 cs_log_printf
2088 (CS_LOG_SETUP,
2089 _(" Large scale Meteo profile info:\n"
2090 " roughness: %12f [m]\n"
2091 " LMO inv: %12f [m^-1]\n"
2092 " ustar0: %12f [m/s]\n"
2093 " uref: %12f [m/s]\n"
2094 " zref: %12f [m]\n"
2095 " angle: %12f [°]\n"
2096 " P sea: %12f [Pa]\n"
2097 " T0: %12f [K]\n"
2098 " Tstar: %12f [K]\n"
2099 " BL height: %12f [m]\n\n"),
2100 cs_glob_atmo_option->meteo_z0,
2101 cs_glob_atmo_option->meteo_dlmo,
2102 cs_glob_atmo_option->meteo_ustar0,
2103 cs_glob_atmo_option->meteo_uref,
2104 cs_glob_atmo_option->meteo_zref,
2105 cs_glob_atmo_option->meteo_angle,
2106 cs_glob_atmo_option->meteo_psea,
2107 cs_glob_atmo_option->meteo_t0,
2108 cs_glob_atmo_option->meteo_tstar,
2109 cs_glob_atmo_option->meteo_zi);
2110 }
2111
2112 }
2113
2114 /*----------------------------------------------------------------------------*/
2115 /*!
2116 * \brief Print the atmospheric chemistry options to setup.log.
2117 */
2118 /*----------------------------------------------------------------------------*/
2119
2120 void
cs_atmo_chemistry_log_setup(void)2121 cs_atmo_chemistry_log_setup(void)
2122 {
2123 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] < 0)
2124 return;
2125
2126 cs_log_printf(CS_LOG_SETUP,
2127 _("\n"
2128 "Atmospheric chemistry options\n"
2129 "-----------------------------\n\n"));
2130
2131 if (cs_glob_atmo_chemistry->model == 0) {
2132
2133 /* No atmospheric chemistry */
2134 cs_log_printf(CS_LOG_SETUP,
2135 _(" No atmospheric chemistry\n\n"));
2136
2137 }
2138 else if ( cs_glob_atmo_chemistry->model == 1
2139 || cs_glob_atmo_chemistry->model == 2
2140 || cs_glob_atmo_chemistry->model == 3) {
2141
2142 /* Pre-defined schemes */
2143 cs_log_printf
2144 (CS_LOG_SETUP,
2145 _(" Atmospheric chemistry activated\n\n"
2146 " Pre-defined scheme %12d\n\n"
2147 " n_species: %18d (Number of species)\n"
2148 " n_reactions: %16d (Number of reactions)\n"
2149 " photolysis: %17s\n"
2150 " frozen_gas_chem: %12s\n\n"),
2151 cs_glob_atmo_chemistry->model,
2152 cs_glob_atmo_chemistry->n_species,
2153 cs_glob_atmo_chemistry->n_reactions,
2154 cs_glob_atmo_chemistry->chemistry_with_photolysis ? "Yes": "No",
2155 cs_glob_atmo_chemistry->frozen_gas_chem ? "Yes": "No");
2156
2157 }
2158 else if (cs_glob_atmo_chemistry->model == 4) {
2159
2160 /* User defined SPACK chemistry */
2161 cs_log_printf
2162 (CS_LOG_SETUP,
2163 _(" Atmospheric chemistry activated\n\n"
2164 " User-defined SPACK scheme\n\n"
2165 " n_species: %18d (Number of species)\n"
2166 " n_reactions: %16d (Number of reactions)\n"
2167 " photolysis: %17s\n"
2168 " frozen_gas_chem: %12s\n"
2169 " Spack file: %17s\n"),
2170 cs_glob_atmo_chemistry->n_species,
2171 cs_glob_atmo_chemistry->n_reactions,
2172 cs_glob_atmo_chemistry->chemistry_with_photolysis ? "Yes": "No",
2173 cs_glob_atmo_chemistry->frozen_gas_chem ? "Yes": "No",
2174 cs_glob_atmo_chemistry->spack_file_name);
2175
2176 }
2177 }
2178
2179 /*----------------------------------------------------------------------------*/
2180 /*!
2181 * \brief Print the atmospheric aerosols options to setup.log.
2182 */
2183 /*----------------------------------------------------------------------------*/
2184
2185 void
cs_atmo_aerosol_log_setup(void)2186 cs_atmo_aerosol_log_setup(void)
2187 {
2188 if (cs_glob_physical_model_flag[CS_ATMOSPHERIC] < 0)
2189 return;
2190
2191 cs_log_printf
2192 (CS_LOG_SETUP,
2193 _("\nAtmospheric aerosols options\n"
2194 "----------------------------\n\n"));
2195
2196 cs_atmo_aerosol_type_t atm_aer_type = cs_glob_atmo_chemistry->aerosol_model;
2197
2198 if (atm_aer_type == CS_ATMO_AEROSOL_OFF)
2199 cs_log_printf(CS_LOG_SETUP,_(" %s\n"),
2200 cs_atmo_aerosol_type_name[atm_aer_type]);
2201
2202 else if (atm_aer_type == CS_ATMO_AEROSOL_SSH) {
2203
2204 /* Atmospheric chemistry with SSH */
2205 cs_log_printf
2206 (CS_LOG_SETUP,
2207 _(" Atmospheric aerosols activated\n\n"
2208 " Global parameters\n\n"
2209 " model: %22s (%s)\n"
2210 " n_layer: %20d (Number of layers inside aerosols)\n"
2211 " n_size: %20d (Number of resolved aerosols)\n"
2212 " init_gas_with_lib: %10s\n"
2213 " init_aero_with_lib: %9s\n"
2214 " namelist: %s\n\n"),
2215 cs_atmo_aerosol_type_enum_name[atm_aer_type],
2216 cs_atmo_aerosol_type_name[atm_aer_type],
2217 cs_glob_atmo_chemistry->n_layer,
2218 cs_glob_atmo_chemistry->n_size,
2219 cs_glob_atmo_chemistry->init_gas_with_lib ? "Yes": "No",
2220 cs_glob_atmo_chemistry->init_aero_with_lib ? "Yes": "No",
2221 cs_glob_atmo_chemistry->aero_file_name);
2222 }
2223 }
2224
2225 /*----------------------------------------------------------------------------*/
2226 /*!
2227 * \brief Deallocate arrays for atmo module
2228 */
2229 /*----------------------------------------------------------------------------*/
2230
2231 void
cs_atmo_finalize(void)2232 cs_atmo_finalize(void)
2233 {
2234 BFT_FREE(_atmo_option.meteo_file_name);
2235 BFT_FREE(_atmo_option.z_temp_met);
2236 BFT_FREE(_atmo_option.time_met);
2237 BFT_FREE(_atmo_option.hyd_p_met);
2238 }
2239
2240 /*----------------------------------------------------------------------------*/
2241
2242 END_C_DECLS
2243