Mass conservation issue for chemical tracers in a continuous run mode about global-workflow

[DWesl]

You mention sophisticated functions to perform vertical interpolation while preserving mass between pressure levels. Could you elaborate on what those functions are?

[zhanglikate]

You mention sophisticated functions to perform vertical interpolation while preserving mass between pressure levels. Could you elaborate on what those functions are?

@rmontuoro is more familiar with the python script vertical interpolation script of https://github.com/NOAA-EMC/global-workflow/blob/develop/scripts/exgfs_aero_init_aerosol.py , who can help to elaborate more details about it.

[DWesl]

It looks like that script passes the processing to another:

global-workflow/scripts/exgfs_aero_init_aerosol.py

Line 52 in b405b7d

merge_script_pattern = "{ush_gfs}/merge_fv3_aerosol_tile.py"

global-workflow/scripts/exgfs_aero_init_aerosol.py

Line 286 in b405b7d

subprocess.run([merge_script, atm_file, tracer_file, core_file, ctrl_file, rest_file, tracer_list_file, temp_file], check=True)

which currently does the aforementioned rescaling here:

global-workflow/ush/merge_fv3_aerosol_tile.py

Lines 70 to 71 in b405b7d

	# read pressure layer thickness from restart file
	delp = rest_file["delp"][0, :]

global-workflow/ush/merge_fv3_aerosol_tile.py

Lines 91 to 93 in b405b7d

	dp = np.zeros(delp.shape)
	for k in range(0, dp.shape[0]):
	dp[k, :, :] = ak[k + 1] - ak[k] + psfc * (bk[k + 1] - bk[k])

global-workflow/ush/merge_fv3_aerosol_tile.py

Line 95 in b405b7d

scale_factor = delp / dp

global-workflow/ush/merge_fv3_aerosol_tile.py

Lines 112 to 113 in b405b7d

	base_file[variable_name][0, :, :] = 0.
	base_file[variable_name][1:, :, :] = scale_factor * variable[0, :, :, :]

It looks like that script rescales the CO₂ mixing ratio to keep the mass per layer constant and sets the layer not available in the restart file to zero.

Given the difference in number of levels between the restart file and the file the model reads in, would setting base_file[variable_name][0, :, :] = base_file[variable_name][1, :, :] after line 113 make sense? If the mass differences printed on lines 105-120 match the mass differences found in the experiment, the zeroed layer is probably not the problem. (EDIT: This is adding a zero on the 1 Pa end of the model, so that in itself would not cause much of a difference, though the lower-CO₂ air aloft being moved surfaceward in the model-internal-to-restart-file transition might still be a concern)

I am reminded that, if the surface pressure changes, it is possible to conserve the CO₂ mass per unit area or mole fraction/mixing ratio values but not both (that is, keeping the mole fraction values at 400ppm will result in the column CO₂ mass changing with the surface pressure).

The issue mentions a proposed solution, described in the accompanying presentation, that achieves better mass conservation than the existing method. I would like more details on this new method.

As an alternative, is there an option to have the model runs be three or five days at a stretch rather than one? That would result in fewer shocks to the system from new meteorological data, with more time as a free run, which doesn't have those problems.

EDIT: Just remembered that runs into problems. Is there an analysis-nudging option?

[yangfanglin]

Daniel, thank you for checking. Looks like the rescaling approach needs further consideration.

I noticed also on slide 3 of Jian He's ppt there is notation Ptop=0. Please note FV3 models cannot have a model top at zero. Ptop should be 20 Pa. If there is a zero pressure in gfs_data.tile[1-6].nc it should have never been used in any model integration. Not sure if this has played any role in interpolation and calculation of mass.