Skip to content
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

Jump to bottom

Disort only tro #508

Merged
merged 4 commits into from
Sep 1, 2022
Merged

Disort only tro #508

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4c872d9
Re-Add removed only_tro-machanism to DisortCalc and DisortCalcWithART…
m-brath Aug 31, 2022
98b1690
Add only_tro-mechanism to DisortCalcIrradiance
m-brath Aug 31, 2022
f16f344
Clean up code
m-brath Aug 31, 2022
6f7df97
Rename DisortCalcClearSky to DisortCalcClearsky
m-brath Aug 31, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion controlfiles/artscomponents/heatingrates/TestSpectralRadianceField.arts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ Copy( ref_field, spectral_radiance_field )

# Set spectral_radiance_field as full DISORT calculation
#
DisortCalcClearSky( nstreams = 8 )
DisortCalcClearsky( nstreams = 8 )
#WriteXML( "binary", spectral_radiance_field, "f2.xml" )
Compare( ref_field, spectral_radiance_field, 2e-18 )
#Copy( ref_field, spectral_radiance_field )
Expand Down
165 changes: 114 additions & 51 deletions src/disort.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -936,37 +936,36 @@ void reduced_1datm(Vector& p,
}

void run_cdisort(Workspace& ws,
// Output
Tensor7& cloudbox_field,
ArrayOfMatrix& disort_aux,
// Input
ConstVectorView f_grid,
ConstVectorView p_grid,
ConstVectorView z_profile,
const Numeric& z_surface,
ConstVectorView t_profile,
ConstMatrixView vmr_profiles,
ConstMatrixView pnd_profiles,
const ArrayOfArrayOfSingleScatteringData& scat_data,
const ArrayOfStar& stars,
const Agenda& propmat_clearsky_agenda,
const Agenda& gas_scattering_agenda,
const ArrayOfIndex& cloudbox_limits,
const Numeric& surface_skin_t,
const Vector& surface_scalar_reflectivity,
ConstVectorView za_grid,
ConstVectorView aa_grid,
ConstVectorView star_rte_los,
const Index& gas_scattering_do,
const Index& stars_do,
const ArrayOfString& disort_aux_vars,
const Numeric& scale_factor,
const Index& nstreams,
const Index& Npfct,
const Index& quiet,
const Index& emission,
const Index& intensity_correction,
const Verbosity& verbosity) {
Tensor7& cloudbox_field,
ArrayOfMatrix& disort_aux,
ConstVectorView f_grid,
ConstVectorView p_grid,
ConstVectorView z_profile,
const Numeric& z_surface,
ConstVectorView t_profile,
ConstMatrixView vmr_profiles,
ConstMatrixView pnd_profiles,
const ArrayOfArrayOfSingleScatteringData& scat_data,
const ArrayOfStar& stars,
const Agenda& propmat_clearsky_agenda,
const Agenda& gas_scattering_agenda,
const ArrayOfIndex& cloudbox_limits,
const Numeric& surface_skin_t,
const Vector& surface_scalar_reflectivity,
ConstVectorView za_grid,
ConstVectorView aa_grid,
ConstVectorView star_rte_los,
const Index& gas_scattering_do,
const Index& stars_do,
const ArrayOfString& disort_aux_vars,
const Numeric& scale_factor,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& quiet,
const Index& emission,
const Index& intensity_correction,
const Verbosity& verbosity) {
// Create an atmosphere starting at z_surface
Vector p, z, t;
Matrix vmr, pnd;
Expand Down Expand Up @@ -1082,21 +1081,53 @@ void run_cdisort(Workspace& ws,
// Transform to mu, starting with negative values
for (Index i = 0; i < ds.numu; i++) ds.umu[i] = -cos(za_grid[i] * PI / 180);

//gas absorption
Matrix ext_bulk_gas(nf, ds.nlyr + 1);
get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);
Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
get_paroptprop(
ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);

// get the angles where to calculate the scattering, which is later used for
//for the calculation of the legendre polynoms
// Get particle bulk properties
Index nang;
Vector pfct_angs;
get_angs(pfct_angs, scat_data, Npfct);
Index nang = pfct_angs.nelem();

Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
Index nf_ssd = scat_data[0][0].f_grid.nelem();
Tensor3 pha_bulk_par(nf_ssd, ds.nlyr + 1, nang);
get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
Tensor3 pha_bulk_par;

if (only_tro && (Npfct < 0 || Npfct > 3)) {
nang = Npfct;
nlinspace(pfct_angs, 0, 180, nang);

pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);

ext_bulk_par = 0.0;
abs_bulk_par = 0.0;
pha_bulk_par = 0.0;

Index iflat = 0;

for (Index iss = 0; iss < scat_data.nelem(); iss++) {
const Index nse = scat_data[iss].nelem();
ext_abs_pfun_from_tro(ext_bulk_par,
abs_bulk_par,
pha_bulk_par,
scat_data[iss],
iss,
pnd(Range(iflat, nse), joker),
cboxlims,
t,
pfct_angs);
iflat += nse;
}
} else {
get_angs(pfct_angs, scat_data, Npfct);
nang = pfct_angs.nelem();

pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);

get_paroptprop(
ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
}

Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);

Expand Down Expand Up @@ -1293,6 +1324,7 @@ void run_cdisort_flux(Workspace& ws,
const Numeric& scale_factor,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& quiet,
const Index& emission,
const Index& intensity_correction,
Expand Down Expand Up @@ -1407,22 +1439,53 @@ void run_cdisort_flux(Workspace& ws,
for (Index i = 0; i <= ds.nlyr; i++) ds.temper[i] = t[ds.nlyr - i];
}

//gas absorption
Matrix ext_bulk_gas(nf, ds.nlyr + 1);
get_gasoptprop(ws, ext_bulk_gas, propmat_clearsky_agenda, t, vmr, p, f_grid);

// Get particle bulk properties
Index nang;
Vector pfct_angs;
Matrix ext_bulk_par(nf, ds.nlyr + 1), abs_bulk_par(nf, ds.nlyr + 1);
get_paroptprop(
ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
Index nf_ssd = scat_data[0][0].f_grid.nelem();
Tensor3 pha_bulk_par;

// get the angles where to calculate the scattering, which is later used for
//for the calculation of the legendre polynoms
Vector pfct_angs;
get_angs(pfct_angs, scat_data, Npfct);
Index nang = pfct_angs.nelem();
if (only_tro && (Npfct < 0 || Npfct > 3)) {
nang = Npfct;
nlinspace(pfct_angs, 0, 180, nang);

pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);

ext_bulk_par = 0.0;
abs_bulk_par = 0.0;
pha_bulk_par = 0.0;

Index iflat = 0;

for (Index iss = 0; iss < scat_data.nelem(); iss++) {
const Index nse = scat_data[iss].nelem();
ext_abs_pfun_from_tro(ext_bulk_par,
abs_bulk_par,
pha_bulk_par,
scat_data[iss],
iss,
pnd(Range(iflat, nse), joker),
cboxlims,
t,
pfct_angs);
iflat += nse;
}
} else {
get_angs(pfct_angs, scat_data, Npfct);
nang = pfct_angs.nelem();

pha_bulk_par.resize(nf_ssd, ds.nlyr + 1, nang);

get_paroptprop(
ext_bulk_par, abs_bulk_par, scat_data, pnd, t, p, cboxlims, f_grid);
get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
}

Index nf_ssd = scat_data[0][0].f_grid.nelem();
Tensor3 pha_bulk_par(nf_ssd, ds.nlyr + 1, nang);
get_parZ(pha_bulk_par, scat_data, pnd, t, pfct_angs, cboxlims);
Tensor3 pfct_bulk_par(nf_ssd, ds.nlyr, nang);
get_pfct(pfct_bulk_par, pha_bulk_par, ext_bulk_par, abs_bulk_par, cboxlims);

Expand Down
6 changes: 4 additions & 2 deletions src/disort.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,7 @@ void get_disortsurf_props( // Output
* @param[in] nstreams Number of quadrature angles (both hemispheres).
* @param[in] Npfct Number of angular grid points to calculate bulk phase
* function.
* @param[in] only_tro Flag indicating that only TRO scattering data is used.
* @param[in] quiet Silence warnings.
* @param[in] emission Enables blackbody emission.
* @param[in] intensity_correction Enables intensity correction (for low nstreams)
Expand All @@ -203,10 +204,8 @@ void get_disortsurf_props( // Output
* @date 2019-09-19, 2021-10-27
*/
void run_cdisort(Workspace& ws,
// Output
Tensor7& cloudbox_field,
ArrayOfMatrix& disort_aux,
// Input
ConstVectorView f_grid,
ConstVectorView p_grid,
ConstVectorView z_profile,
Expand All @@ -230,6 +229,7 @@ void run_cdisort(Workspace& ws,
const Numeric& scale_factor,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& quiet,
const Index& emission,
const Index& intensity_correction,
Expand Down Expand Up @@ -274,6 +274,7 @@ void run_cdisort(Workspace& ws,
* @param[in] nstreams Number of quadrature angles (both hemispheres).
* @param[in] Npfct Number of angular grid points to calculate bulk phase
* function.
* @param[in] only_tro Flag indicating that only TRO scattering data is used.
* @param[in] quiet Silence warnings.
* @param[in] emission Enables blackbody emission.
* @param[in] intensity_correction Enables intensity correction (for low nstreams)
Expand Down Expand Up @@ -306,6 +307,7 @@ void run_cdisort_flux(Workspace& ws,
const Numeric& scale_factor,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& quiet,
const Index& emission,
const Index& intensity_correction,
Expand Down
28 changes: 21 additions & 7 deletions src/m_disort.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,7 @@ void DisortCalc(Workspace& ws,
const ArrayOfString& disort_aux_vars,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& cdisort_quiet,
const Index& emission,
const Index& intensity_correction,
Expand Down Expand Up @@ -147,7 +148,9 @@ void DisortCalc(Workspace& ws,
if (star_rte_los[0] >= 90) {
star_on = 0;

//TODO: Add warning message that star is switched off because it is below horizon
CREATE_OUT0;
out0 << "Star is below the horizon\n";
out0 << "Star is ignored.\n";
}

init_ifield(cloudbox_field,
Expand Down Expand Up @@ -180,6 +183,10 @@ void DisortCalc(Workspace& ws,


} else {
CREATE_OUT3;
out3 << "Disort calculation encountered aa_grid size larger than 1 in a case when it\n";
out3 << "does not use aa_grid. Calculations are performed as if there is no aa_grid.\n";

init_ifield(cloudbox_field,
f_grid,
cloudbox_limits,
Expand Down Expand Up @@ -220,6 +227,7 @@ void DisortCalc(Workspace& ws,
scale_factor,
nstreams,
Npfct,
only_tro,
cdisort_quiet,
emission,
intensity_correction,
Expand Down Expand Up @@ -263,6 +271,7 @@ void DisortCalcWithARTSSurface(Workspace& ws,
const ArrayOfString& disort_aux_vars,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& cdisort_quiet,
const Index& emission,
const Index& intensity_correction,
Expand Down Expand Up @@ -332,7 +341,9 @@ void DisortCalcWithARTSSurface(Workspace& ws,
if (star_rte_los[0] >= 90) {
star_on = 0;

//TODO: Add warning message that star is switched off because it is below horizon
CREATE_OUT0;
out0 << "Star is below the horizon\n";
out0 << "Star is ignored.\n";
}

init_ifield(cloudbox_field,
Expand Down Expand Up @@ -425,6 +436,7 @@ void DisortCalcWithARTSSurface(Workspace& ws,
scale_factor,
nstreams,
Npfct,
only_tro,
cdisort_quiet,
emission,
intensity_correction,
Expand All @@ -433,7 +445,7 @@ void DisortCalcWithARTSSurface(Workspace& ws,


/* Workspace method: Doxygen documentation will be auto-generated */
void DisortCalcClearSky(Workspace& ws,
void DisortCalcClearsky(Workspace& ws,
// WS Output:
Tensor7& spectral_radiance_field,
ArrayOfMatrix& disort_aux,
Expand Down Expand Up @@ -540,6 +552,7 @@ void DisortCalcClearSky(Workspace& ws,
nstreams,
181,
cdisort_quiet,
0,
emission,
intensity_correction,
verbosity);
Expand Down Expand Up @@ -578,17 +591,15 @@ void DisortCalcIrradiance(Workspace& ws,
const ArrayOfString& disort_aux_vars,
const Index& nstreams,
const Index& Npfct,
const Index& only_tro,
const Index& cdisort_quiet,
const Index& emission,
const Index& intensity_correction,
const Verbosity& verbosity) {
// Don't do anything if there's no cloudbox defined.

// Set cloudbox to cover complete atmosphere
Index cloudbox_on;
ArrayOfIndex cloudbox_limits;
const Index cloudbox_checked = 1;
//
cloudboxSetFullAtm(cloudbox_on,
cloudbox_limits,
atmosphere_dim,
Expand Down Expand Up @@ -656,7 +667,9 @@ void DisortCalcIrradiance(Workspace& ws,
if (star_rte_los[0] >= 90) {
star_on = 0;

//TODO: Add warning message that star is switched off because it is below horizon
CREATE_OUT0;
out0 << "Star is below the horizon\n";
out0 << "Star is ignored.\n";
}

//get the cloudbox top distance to earth center.
Expand Down Expand Up @@ -711,6 +724,7 @@ void DisortCalcIrradiance(Workspace& ws,
scale_factor,
nstreams,
Npfct,
only_tro,
cdisort_quiet,
emission,
intensity_correction,
Expand Down
Loading