Skip to content

Aniso smg #172

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pdziekan wants to merge 17 commits into
base: master
Choose a base branch
from
Open

Aniso smg #172

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
be2d05d
anisotropic SMG with mixing lengths dependant on resolution along the…
pdziekan Feb 14, 2024
58d5af3
k_ma limiters
pdziekan Feb 16, 2024
dc6a24a
Merge branch 'blk_trapez_fixes' into aniso_smg
pdziekan Feb 16, 2024
d3d61c0
Merge branch 'master' into aniso_smg
pdziekan Feb 21, 2024
dc6d008
Merge branch 'master' into aniso_smg
pdziekan Nov 21, 2024
f092d8b
central gradient in subsidence instead of forward to fix top boundary…
Dec 2, 2024
b47512f
fix boost phoenix linkage for boost 1.81-1.83
Dec 20, 2024
5a7e8cb
Merge branch 'master' into aniso_smg
May 14, 2025
98b32c0
comment
pdziekan Dec 9, 2025
d8b0144
Merge branch 'master' of github.com:igfuw/UWLCM into aniso_smg
pdziekan Jan 8, 2026
37a5e0d
run_hlpr: sgs as string option accepting iles,smg,smagni
pdziekan Jan 8, 2026
65a2b48
smganiso WIP
pdziekan Jan 12, 2026
d0f5882
slvr_sgs: smg_common class
pdziekan Jan 12, 2026
7f7c425
Merge branch 'master' of github.com:igfuw/UWLCM into aniso_smg
pdziekan Jan 12, 2026
1bb554f
smgani solver
pdziekan Jan 19, 2026
488abcc
smgani unit test
pdziekan Jan 19, 2026
8e8b257
GA test smgani fixes
pdziekan Jan 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
The table of contents is too big for display.
Diff view
Diff view
  •  
  •  
  •  
The diff you're trying to view is too large. We only load the first 3000 changed files.
11 changes: 7 additions & 4 deletions .github/workflows/test_uwlcm_hlpr.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,18 +74,21 @@ jobs:
matrix:
build_type: ["RelWithDebInfo", "Debug"]
mpi: ["none"]
sgs: ["iles", "smg"]
sgs: ["iles", "smg", "smgani"]
disable_cuda: ["true"]
include:
- mpi: "none"
tag: "uwlcm_ubuntu_24_04_cuda_12_9_0"
cxx: "g++"
- sgs: "iles"
cmake_flags: "-DUWLCM_DISABLE='SGS'"
cmake_flags: "-DUWLCM_DISABLE='SMG;SMGANI'"
ctests: "api_test_iles"
- sgs: "smg"
cmake_flags: "-DUWLCM_DISABLE='ILES'"
cmake_flags: "-DUWLCM_DISABLE='ILES;SMGANI'"
ctests: "api_test_smg"
- sgs: "smgani"
cmake_flags: "-DUWLCM_DISABLE='ILES;SMG'"
ctests: "api_test_smgani"
- build_type: "RelWithDebInfo"
tar_name: "RelDbgPort"
- build_type: "Debug"
Expand Down Expand Up @@ -231,7 +234,7 @@ jobs:
- name: configure and make UWLCM
working-directory: ${{ github.workspace }}
run: VERBOSE=1 OMP_NUM_THREADS=4 singularity exec $SI bash -c "
cmake -B build -DCMAKE_BUILD_TYPE=${{matrix.build_type}} -DCMAKE_CXX_COMPILER=${{matrix.cxx}} -Dlibcloudph++_DIR=${{ github.workspace }}/installed/share/libcloudph++ -Dlibmpdata++_DIR=${{ github.workspace }}/installed/share/libmpdata++ -DUWLCM_DISABLE='3D_LGRNGN;3D_BLK_1M;3D_BLK_2M;PIGGYBACKER;SGS' -DMPDATA_OPTS='ABS' -DUWLCM_TIMING=1 -DCMAKE_INSTALL_PREFIX=${{ github.workspace }}/installed &&
cmake -B build -DCMAKE_BUILD_TYPE=${{matrix.build_type}} -DCMAKE_CXX_COMPILER=${{matrix.cxx}} -Dlibcloudph++_DIR=${{ github.workspace }}/installed/share/libcloudph++ -Dlibmpdata++_DIR=${{ github.workspace }}/installed/share/libmpdata++ -DUWLCM_DISABLE='3D_LGRNGN;3D_BLK_1M;3D_BLK_2M;PIGGYBACKER;SMG;SMGANI' -DMPDATA_OPTS='ABS' -DUWLCM_TIMING=1 -DCMAKE_INSTALL_PREFIX=${{ github.workspace }}/installed &&
cmake --build build --config ${{matrix.build_type}} -j2
"

Expand Down
10 changes: 6 additions & 4 deletions src/cases/BOMEX03.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -237,7 +237,7 @@ namespace cases
// like in Wojtek's BabyEulag
// alse set w_LS and hgt_fctrs
// TODO: same in DYCOMS (and others?), move to a common function
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
using libcloudphxx::common::moist_air::R_d_over_c_pd;
using libcloudphxx::common::moist_air::c_pd;
Expand All @@ -246,9 +246,10 @@ namespace cases
using libcloudphxx::common::theta_std::p_1000;

blitz::firstIndex k;
const int nz = nps[n_dims - 1];
real_t dz = (this->Z / si::metres) / (nz-1);

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
parent_t::env_prof(profs, nz);
parent_t::ref_prof(profs, nz);

Expand Down Expand Up @@ -414,11 +415,12 @@ namespace cases
concurr.vab_relaxed_state(1) = concurr.advectee(ix::v);
}

void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[3], const user_params_t &user_params)
{
parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
// geostrophic wind used in Coriolis force
blitz::firstIndex k;
const int nz = nps[2];
real_t dz = (this->Z / si::metres) / (nz-1);
profs.geostr[0] = -10 + 1.8e-3 * k * dz;
profs.geostr[1] = 0;
Expand Down
65 changes: 59 additions & 6 deletions src/cases/CasesCommon.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,9 @@ namespace cases
template<class case_ct_params_t, int n_dims>
class CasesCommon
{
static constexpr real_t karman_c = 0.41; // von Karman constant
static constexpr real_t smg_c = 0.165; // Smagorinsky constant

public:

using ix = typename case_ct_params_t::ix;
Expand Down Expand Up @@ -102,8 +105,9 @@ namespace cases
typename std::enable_if<enable_sgs>::type* = 0)
{
params.c_m = 0.0856;
params.smg_c = 0.165;
params.smg_c = smg_c;
params.prandtl_num = 0.42;
params.karman_c = karman_c;
params.cdrag = 0;
params.fricvelsq = 0;
}
Expand All @@ -113,23 +117,72 @@ namespace cases

virtual void setopts(rt_params_t &params, const int nps[], const user_params_t &user_params) {assert(false);};
virtual void intcond(concurr_any_t &concurr, arr_1D_t &rhod, arr_1D_t &th_e, arr_1D_t &rv_e, arr_1D_t &rl_e, arr_1D_t &p_e, int rng_seed) =0;
virtual void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)

// virtual void set_profs(detail::profiles_t &profs, int nx, int nz, const user_params_t &user_params)
virtual void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
real_t dz = (Z / si::metres) / (nz-1);
static_assert(n_dims == 2 || n_dims == 3, "set_profs: n_dims must be 2 or 3");
// TODO: get these constants from user_params (or rt_params_t)
// const real_t smg_c = 0.165; // Smagornsky constant
// const real_t karman_c = 0.41; // von Karman constant

const int nz = nps[n_dims-1];

// grid spacing
real_t dx, dy, dz;
if constexpr (n_dims == 2)
{
dx = (X / si::metres) / (nps[0]-1);
dy = dx; // used in horizontal mixing length...
dz = (Z / si::metres) / (nps[1]-1);
}
else if (n_dims == 3)
{
dx = (X / si::metres) / (nps[0]-1);
dy = (Y / si::metres) / (nps[1]-1);
dz = (Z / si::metres) / (nps[2]-1);
}

// set SGS mixing length
{
using blitz::pow2;
blitz::firstIndex k;

real_t sgs_delta;
real_t sgs_delta; // size of a grid cell used as characteristic length scale in SGS models
if (user_params.sgs_delta > 0)
{
sgs_delta = user_params.sgs_delta;
}
else
{
sgs_delta = dz;
sgs_delta = dz; // default
}
profs.mix_len = min(max(k, 1) * dz * 0.845, sgs_delta);

// size of SGS eddies, used in SD SGS models (note that they are isotropic!)
profs.sgs_eddy_size = min(max(k, 0.5) * dz * 0.845, sgs_delta);

// isotropic mixing length - used in isotropic SMG and in SD SGS models
// like in Smagorinsky 1963, but with default characteristic length scale (sgs_delta) equal to dz (not dx*dy*dz)^(1/3))
// similar approach used e.g. in the SAM model
// also there's a reduction near ground
// NOTE: it's not multiplied by the smagorinsky constant here (that is done in slvr_sgs), because mix_len_iso is used also in SD SGS (opts_lgrngn.hpp) where no such constant is used
// profs.mix_len_iso = min(max(k, 0.5) * dz * 0.845, sgs_delta);
// profs.mix_len_iso_sq = pow2(smg_c * min(max(k, 0.5) * dz * 0.845, sgs_delta));
Comment on lines +161 to +170
Copy link

Copilot AI Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The magic number 0.845 appears multiple times (lines 162, 171, 180, 181) without explanation. While there's a comment on line 167 mentioning reduction near ground, it would be helpful to document what this constant represents (appears to be related to von Kármán constant or mixing length reduction factor).

Suggested change
// size of SGS eddies, used in SD SGS models (note that they are isotropic!)
profs.sgs_eddy_size = min(max(k, 0.5) * dz * 0.845, sgs_delta);
// isotropic mixing length - used in isotropic SMG and in SD SGS models
// like in Smagorinsky 1963, but with default characteristic length scale (sgs_delta) equal to dz (not dx*dy*dz)^(1/3))
// similar approach used e.g. in the SAM model
// also there's a reduction near ground
// NOTE: it's not multiplied by the smagorinsky constant here (that is done in slvr_sgs), because mix_len_iso is used also in SD SGS (opts_lgrngn.hpp) where no such constant is used
// profs.mix_len_iso = min(max(k, 0.5) * dz * 0.845, sgs_delta);
// profs.mix_len_iso_sq = pow2(smg_c * min(max(k, 0.5) * dz * 0.845, sgs_delta));
// empirical factor reducing SGS eddy size / mixing length near the ground
const real_t sgs_near_ground_reduction = real_t(0.845);
// size of SGS eddies, used in SD SGS models (note that they are isotropic!)
profs.sgs_eddy_size = min(max(k, 0.5) * dz * sgs_near_ground_reduction, sgs_delta);
// isotropic mixing length - used in isotropic SMG and in SD SGS models
// like in Smagorinsky 1963, but with default characteristic length scale (sgs_delta) equal to dz (not dx*dy*dz)^(1/3))
// similar approach used e.g. in the SAM model
// also there's a reduction near ground (see sgs_near_ground_reduction above)
// NOTE: it's not multiplied by the smagorinsky constant here (that is done in slvr_sgs), because mix_len_iso is used also in SD SGS (opts_lgrngn.hpp) where no such constant is used
// profs.mix_len_iso = min(max(k, 0.5) * dz * sgs_near_ground_reduction, sgs_delta);
// profs.mix_len_iso_sq = pow2(smg_c * min(max(k, 0.5) * dz * sgs_near_ground_reduction, sgs_delta));

Copilot uses AI. Check for mistakes.
profs.mix_len_iso_sq = pow2(smg_c * profs.sgs_eddy_size);

// anisotropic mixing length - used in anisotropic SMG
// based on Mason & Thomson 1992, but anisotropic, like in Simon & Chow 2021 (https://doi.org/10.1007/s10546-021-00642-0)
// but we do not use grid size in the normal direction, only along (i.e. dx, dy for horizontal and dz for vertical)
// TODO: to be consistent, setting user_params.sgs_delta or turning on SD SGS models should not be allowed with anisotropic SMG (?)
// NOTE: it's not really mixing length, but grid size used to calculate mixing length; mixing length calculation is in slvr_sgs.hpp;
// here we just set the characteristic grid size along each direction to simplify potential future option to override it (as with isotropic sgs_delta)
real_t dxy = pow(dx * dy, real_t(1./2));
profs.mix_len_hori_sq = pow2(real_t(1) / (real_t(1) / pow2(karman_c * max(k, 1) * dz) + real_t(1) / pow2(smg_c * dxy)));
profs.mix_len_vert_sq = pow2(real_t(1) / (real_t(1) / pow2(karman_c * max(k, 1) * dz) + real_t(1) / pow2(smg_c * dz)));
Comment on lines +180 to +181
Copy link

Copilot AI Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The mixing length formula is complex and spans one line. Consider breaking it into intermediate variables or adding a comment explaining the formula structure (appears to be harmonic mean of two mixing lengths).

Suggested change
profs.mix_len_hori_sq = pow2(real_t(1) / (real_t(1) / pow2(karman_c * max(k, 1) * dz) + real_t(1) / pow2(smg_c * dxy)));
profs.mix_len_vert_sq = pow2(real_t(1) / (real_t(1) / pow2(karman_c * max(k, 1) * dz) + real_t(1) / pow2(smg_c * dz)));
// combine wall-based and grid-based length scales using a harmonic mean of their squared values
real_t wall_mix_len = karman_c * max(k, 1) * dz;
real_t hori_grid_mix_len = smg_c * dxy;
real_t vert_grid_mix_len = smg_c * dz;
real_t inv_wall_len_sq = real_t(1) / pow2(wall_mix_len);
real_t inv_hori_grid_len_sq = real_t(1) / pow2(hori_grid_mix_len);
real_t inv_vert_grid_len_sq = real_t(1) / pow2(vert_grid_mix_len);
real_t hori_mix_len = real_t(1) / (inv_wall_len_sq + inv_hori_grid_len_sq);
real_t vert_mix_len = real_t(1) / (inv_wall_len_sq + inv_vert_grid_len_sq);
profs.mix_len_hori_sq = pow2(hori_mix_len);
profs.mix_len_vert_sq = pow2(vert_mix_len);

Copilot uses AI. Check for mistakes.
// real_t(1) / (real_t(1) / pow(this->smg_c * params.di, 2) + (*this->params.mix_len)(this->vert_idx))
// profs.mix_len_vert_sq = dz; // real_t(1) / pow2(karman_c * max(k, 1) * dz);
// add dx dy terms; also add sgs_delta, from user params...
// TODO: save these, use them in slvr_sgs, ...
Copy link

Copilot AI Jan 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These commented-out lines and TODOs should be removed or implemented. They appear to be development notes that should not remain in production code.

Suggested change
// TODO: save these, use them in slvr_sgs, ...

Copilot uses AI. Check for mistakes.
}

// set profile of vertical distribution of surface fluxes
Expand Down
6 changes: 4 additions & 2 deletions src/cases/CumulusCongestus_icmw20.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,15 +37,17 @@ namespace cases

// calculate the initial environmental theta and rv profiles
// alse set w_LS and hgt_fctrs
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
using libcloudphxx::common::moist_air::R_d_over_c_pd;
using libcloudphxx::common::moist_air::c_pd;
using libcloudphxx::common::moist_air::R_d;
using libcloudphxx::common::const_cp::l_tri;
using libcloudphxx::common::theta_std::p_1000;

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);

const int nz = nps[n_dims - 1];

// read the soundings
// containers for soundings
Expand Down
5 changes: 3 additions & 2 deletions src/cases/CumulusCongestus_icmw24.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,10 +144,11 @@ namespace cases

// calculate the initial environmental theta and rv profiles
// alse set w_LS and hgt_fctrs
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);

const int nz = nps[n_dims - 1];
this->env_prof(profs, nz);
this->ref_prof(profs, nz);

Expand Down
11 changes: 6 additions & 5 deletions src/cases/DYCOMS.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -229,7 +229,7 @@ namespace cases
// like in Wojtek's BabyEulag
// alse set w_LS and hgt_fctrs
// TODO: move hgt_fctrs from cases to main code
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params) override
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params) override
{
using libcloudphxx::common::moist_air::R_d_over_c_pd;
using libcloudphxx::common::moist_air::c_pd;
Expand All @@ -238,9 +238,10 @@ namespace cases
using libcloudphxx::common::theta_std::p_1000;

blitz::firstIndex k;
const int nz = nps[n_dims - 1];
real_t dz = (this->Z / si::metres) / (nz-1);

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
parent_t::env_prof(profs, nz);
parent_t::ref_prof(profs, nz);

Expand Down Expand Up @@ -403,12 +404,12 @@ namespace cases
concurr.vab_relaxed_state(1) = concurr.advectee(ix::v);
}

void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[3], const user_params_t &user_params)
{
parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
// geostrophic wind equal to the initial velocity profile
blitz::firstIndex k;
real_t dz = (this->Z / si::metres) / (nz-1);
real_t dz = (this->Z / si::metres) / (nps[2]-1);
profs.geostr[0] = this->u(k * dz);
profs.geostr[1] = v(k * dz);
}
Expand Down
5 changes: 3 additions & 2 deletions src/cases/DryPBL.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,12 +224,13 @@ namespace cases
}
}

void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
blitz::firstIndex k;
const int nz = nps[n_dims - 1];
real_t dz = (this->Z / si::metres) / (nz-1);

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
parent_t::env_prof(profs, nz);
parent_t::ref_prof(profs, nz);

Expand Down
4 changes: 2 additions & 2 deletions src/cases/DryThermalGMD2015.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,14 +83,14 @@ namespace cases
// calculate the initial environmental theta and rv profiles
// alse set w_LS and hgt_fctrs
// like in Wojtek's BabyEulag
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
using libcloudphxx::common::theta_std::p_1000;
using libcloudphxx::common::moist_air::R_d;
using libcloudphxx::common::moist_air::c_pd;
using libcloudphxx::common::moist_air::R_d_over_c_pd;

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);

profs.rhod = 1;
profs.th_e = 300;
Expand Down
5 changes: 3 additions & 2 deletions src/cases/MoistThermalGrabowskiClark99.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -192,7 +192,7 @@ namespace cases
public:
// calculate the initial environmental theta and rv profiles as Wojtek does it
// i.e. for stable virtual standard potential temperature
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
// pre_ref - total pressure
// th_e - dry potential temp
// th_ref - dry potential temp refrence profsile
Expand All @@ -206,8 +206,9 @@ namespace cases
using libcloudphxx::common::theta_std::p_1000;
using setup::real_t;

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);

const int nz = nps[n_dims - 1];
real_t dz = (this->Z / si::metres) / (nz-1);
blitz::firstIndex k;
// temperature and total pressure profiles
Expand Down
11 changes: 6 additions & 5 deletions src/cases/RICO11.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,7 +284,7 @@ namespace cases
// like in Wojtek's BabyEulag
// alse set w_LS and hgt_fctrs
// TODO: same in DYCOMS (and others?), move to a common function
void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[n_dims], const user_params_t &user_params)
{
using libcloudphxx::common::moist_air::R_d_over_c_pd;
using libcloudphxx::common::moist_air::c_pd;
Expand All @@ -293,9 +293,10 @@ namespace cases
using libcloudphxx::common::theta_std::p_1000;

blitz::firstIndex k;
const int nz = nps[n_dims - 1];
real_t dz = (this->Z / si::metres) / (nz-1);

parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
parent_t::env_prof(profs, nz);
parent_t::ref_prof(profs, nz);

Expand Down Expand Up @@ -456,12 +457,12 @@ namespace cases
concurr.vab_relaxed_state(1) = concurr.advectee(ix::v);
}

void set_profs(detail::profiles_t &profs, int nz, const user_params_t &user_params)
void set_profs(detail::profiles_t &profs, const int nps[3], const user_params_t &user_params)
{
parent_t::set_profs(profs, nz, user_params);
parent_t::set_profs(profs, nps, user_params);
// geostrophic wind equal to the initial velocity profile, doesnt include mean, because its only used in coriolis = u-geostr
blitz::firstIndex k;
real_t dz = (this->Z / si::metres) / (nz-1);
real_t dz = (this->Z / si::metres) / (nps[2]-1);
profs.geostr[0] = this->u(k * dz);
profs.geostr[1] = v(k * dz);
}
Expand Down
11 changes: 7 additions & 4 deletions src/detail/profiles.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,13 @@ namespace detail
// TODO: store profiles in mem (add an alloc 1d function to libmpdata?) and not in params
struct profiles_t
{
setup::arr_1D_t th_e, p_e, rv_e, rl_e, th_ref, rhod, w_LS, hgt_fctr, th_LS, rv_LS, mix_len, relax_th_rv_coeff;
setup::arr_1D_t th_e, p_e, rv_e, rl_e, th_ref, rhod, w_LS, hgt_fctr, th_LS, rv_LS, mix_len_iso_sq, mix_len_hori_sq, mix_len_vert_sq, sgs_eddy_size, relax_th_rv_coeff;
std::array<setup::arr_1D_t, 2> geostr;

profiles_t(int nz) :
// rhod needs to be bigger, cause it divides vertical courant number
// TODO: should have a halo both up and down, not only up like now; then it should be interpolated in courant calculation
th_e(nz), p_e(nz), rv_e(nz), rl_e(nz), th_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr(nz), th_LS(nz), rv_LS(nz), mix_len(nz), relax_th_rv_coeff(nz)
th_e(nz), p_e(nz), rv_e(nz), rl_e(nz), th_ref(nz), rhod(nz+1), w_LS(nz), hgt_fctr(nz), th_LS(nz), rv_LS(nz), mix_len_iso_sq(nz), mix_len_hori_sq(nz), mix_len_vert_sq(nz), sgs_eddy_size(nz), relax_th_rv_coeff(nz)
{
geostr[0].resize(nz);
geostr[1].resize(nz);
Expand All @@ -32,7 +32,7 @@ namespace detail

struct profile_ptrs_t
{
setup::arr_1D_t *th_e, *p_e, *rv_e, *rl_e, *th_ref, *rhod, *w_LS, *hgt_fctr, *geostr[2], *th_LS, *rv_LS, *mix_len, *relax_th_rv_coeff;
setup::arr_1D_t *th_e, *p_e, *rv_e, *rl_e, *th_ref, *rhod, *w_LS, *hgt_fctr, *geostr[2], *th_LS, *rv_LS, *mix_len_iso_sq, *mix_len_hori_sq, *mix_len_vert_sq, *sgs_eddy_size, *relax_th_rv_coeff;
};

// copy external profiles into rt_parameters
Expand All @@ -53,7 +53,10 @@ namespace detail
{p.rv_LS , profs.rv_LS },
{p.geostr[0] , profs.geostr[0] },
{p.geostr[1] , profs.geostr[1] },
{p.mix_len , profs.mix_len },
{p.mix_len_iso_sq , profs.mix_len_iso_sq },
{p.mix_len_hori_sq , profs.mix_len_hori_sq },
{p.mix_len_vert_sq , profs.mix_len_vert_sq },
{p.sgs_eddy_size , profs.sgs_eddy_size },
{p.relax_th_rv_coeff , profs.relax_th_rv_coeff }
};

Expand Down
2 changes: 1 addition & 1 deletion src/opts/opts_lgrngn.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ void setopts_micro(
neg_w_LS *= -1.; // libcloudphxx defines w_LS>0 for downward direction
std::vector<setup::real_t> vneg_w_LS(neg_w_LS.begin(), neg_w_LS.end());
rt_params.cloudph_opts_init.w_LS = vneg_w_LS;
rt_params.cloudph_opts_init.SGS_mix_len = std::vector<setup::real_t>(rt_params.mix_len->begin(), rt_params.mix_len->end());
rt_params.cloudph_opts_init.SGS_mix_len = std::vector<setup::real_t>(rt_params.sgs_eddy_size->begin(), rt_params.sgs_eddy_size->end());
rt_params.cloudph_opts_init.aerosol_independent_of_rhod = rt_params.aerosol_independent_of_rhod;
rt_params.cloudph_opts_init.aerosol_conc_factor = rt_params.aerosol_conc_factor;

Expand Down
Loading
Loading