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.

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

Add variable rescaling for aerosol species #1339

Draft
wants to merge 10 commits into
base: develop
Choose a base branch
from
21 changes: 20 additions & 1 deletion parm/aero/berror/aero_diagb.yaml.j2
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,26 @@ variables:
- mass_fraction_of_sea_salt003_in_air
- mass_fraction_of_sea_salt004_in_air

rescale: {{ aero_diagb_rescale }} # rescales the filtered std. dev. by "rescale"
global rescale:
geometry:
fms initialization:
namelist filename: ./fv3jedi/fmsmpp.nml
field table filename: ./fv3jedi/field_table
akbk: ./fv3jedi/akbk.nc4
layout:
- {{ layout_x }}
- {{ layout_y }}
nxp: {{ npx_rescale }}
npy: {{ npy_rescale }}
npz: {{ npz_ges }}
field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml
rescale stddev:
filetype: fms restart
skip coupler file: true
dateapath: ./stddev
filename_trcr: rescale.fv_tracer.res.nc
filename_cplr: rescale.coupler.res

number of halo points: {{ aero_diagb_n_halo }}
number of neighbors: {{ aero_diagb_n_neighbors }}
simple smoothing:
Expand Down
42 changes: 37 additions & 5 deletions utils/chem/chem_diagb.h
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,11 @@ namespace gdasapp {
oops::Log::info() << "Background:" << std::endl;
oops::Log::info() << xb << std::endl;

/// Read the global rescale
//--------------------------
// oops::Log::info() << "=================== read global rescale" << std::endl;
// fv3jedi::

/// Create the mesh connectivity (Copy/paste of Francois's stuff)
// --------------------------------------------------------------
// Build edges, then connections between nodes and edges
Expand Down Expand Up @@ -230,13 +235,39 @@ namespace gdasapp {
}
}

// Rescale
if (fullConfig.has("rescale")) {
double rescale;
fullConfig.get("rescale", rescale);
util::multiplyFieldSet(bkgErrFs, rescale);
Copy link
Contributor

Choose a reason for hiding this comment

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

What about using multiplyFieldSets? I think this would require creating a 3D array of values in a new field set, but that would then set us up for the ability to scale by height/latitude/location at a later date. Thoughts?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I can look at it and see how that would work. It took me some time to make something that could run....

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I did try multiplyFieldSets, but had trouble getting it to compile.

Copy link
Contributor

Choose a reason for hiding this comment

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

I wonder if we just read in another 'state' that is the scaling factors and be proactive on the 3D front?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

So have a file with the scaling factors, and read it in from there?

Copy link
Contributor

Choose a reason for hiding this comment

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

I think so, that way we can easily scale by lat,lon,height and species

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Oh, I think I understand. A full background state file, filled with rescaling values for each aerosol species and at all locations. I'd like to try just 14 values first and see if we can do better.

Copy link
Contributor

Choose a reason for hiding this comment

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

yeah they can be constant in x,y,z to start

// Rescale
if (fullConfig.has("global rescale")) {
const eckit::LocalConfiguration GlobalRescaleConfig(fullConfig, "global rescale");
// fv3jedi::State global_rescale(rescaleGeom, chemVars,cycleDate);
// fv3jedi::State global_rescale_interp(geom, chemVars, cycleDate);
const eckit::LocalConfiguration GlobalRescaleGeomConfig(GlobalRescaleConfig,"geometry");
const fv3jedi::Geometry GlobalRescaleGeom(GlobalRescaleGeomConfig, this-> getComm());
fv3jedi::Increment global_rescale(GlobalRescaleGeom, chemVars, cycleDate);
global_rescale.zero();
const eckit::LocalConfiguration GlobalRescaleStdConfig(GlobalRescaleConfig,"rescale stddev");
// Get the 'datapath' and 'filename_trcr' from the YAML configuration
std::string datapath, filename_trcr;
GlobalRescaleStdConfig.get("datapath", datapath);
GlobalRescaleStdConfig.get("filename_trcr", filename_trcr);

// Combine the path and filename to get the full file path
std::string fullPath = datapath + filename_trcr;

// Print out the full path to verify
std::cout << "Attempting to read the global rescale file from: " << fullPath << std::endl;

global_rescale.read(GlobalRescaleStdConfig);
// interpolate to background resolution
fv3jedi::Increment global_rescale_interp(geom, global_rescale);
atlas::FieldSet xrsFs;
global_rescale_interp.toFieldSet(xrsFs);
oops::Log::info() << "global rescaling coefficients:" << std::endl;
oops::Log::info() << xrsFs << std::endl;
util::multiplyFieldSets(bkgErrFs, xrsFs);
}



bkgErr.fromFieldSet(bkgErrFs);

// Hybrid B option
Expand Down Expand Up @@ -274,6 +305,7 @@ namespace gdasapp {
double rescale_staticb;
ClimBConfig.get("staticb rescaling factor", rescale_staticb);


// Combine diagb and climatological background errors
fv3jedi::Increment stddev_hybrid(geom, chemVars, cycleDate);
stddev_hybrid.zero();
Expand Down
Loading