Modify the surface air pressure in RCE model based on climt

[MengZuo]

• CliMT version: climt (v. 0.15.3)
• Python version: 3.10.12
• Operating System:Linux 2.6.32-754.10.1.el6.x86_64

Description

I have recently been learning how to use the RCE model and encountered some issues when trying to modify the surface air pressure.

My purpose is to set the surface air pressure to 700hPa and analyze the elevated heating effect, just like Hu and Boos (2017) did, they used the CliMT (Caballero 2012).

Hu, S., and W. R. Boos, 2017: The Physics of Orographic Elevated Heating in Radiative–Convective Equilibrium. J. Atmos. Sci., 74, 2949–2965, https://doi.org/10.1175/JAS-D-16-0312.1.

My current approach is to adjust the surface air pressure based on gmd_radiative_convective.py. However, I encounter “Segmentation fault (core dumped)” issues when running it.

Do you have suggestions on how to modify the surface air pressure in this model?

I added the following line in gmd_radiative_convective.py:
state['surface_air_pressure'].values[:] = 70000.0

errors are shown here:
19 [[295.72299787]]
[[4201.48460143]]
[[18479.28398476]]
Segmentation fault (core dumped)

[JoyMonteiro]

Thanks for the report! Could you confirm that you are actually using version 0.15.3 ? There have been many updates to climt after that.

[MengZuo]

Hi Joy, actually I'm not sure which version I'm using, but I saw it was updated on February 1, 2023. I installed it in July using pip install, so I believe it should be the latest version.

[MengZuo]

Hi Joy, do you have any ideas about modifying the surface air pressure in RCE model? Thanks!

[JoyMonteiro]

Sorry for the delayed reply, have been travelling.

I think the issue is due to the boundary layer component, which is quite simple.

if you see the initialization function

climt/climt/_components/simple_physics/component.py

Line 97 in e6e257d

top_of_boundary_layer=85000.0,

you will see that the top of the boundary layer is set to 850 hPa, and setting the surface pressure to below this might be causing the issue that you observe. My suggestion is to set this parameter to maybe 700 hPa (you can try a few options depending on your observations).

Also, to ensure the boundary layer is not too deep, maybe you could also play around with the boundary layer influence height

climt/climt/_components/simple_physics/component.py

Line 98 in e6e257d

boundary_layer_influence_height=20000.0,

Maybe reduce this by half to ensure the boundary layer is not too deep.

Also, you should be setting the surface pressure while creating the model grid. Since climt uses a hybrid sigma grid, the surface pressure is an important part of the grid itself and climt must know about it while creating model state. So, if you make the following change

grid_state = get_grid(p_surf_in_Pa=70000.0)

state = get_default_state([simple_physics, convection,
                           radiation_lw, radiation_sw, slab], grid_state=grid_state)

If I make these two changes, the model runs fine. Please let me know if this solves your problem!

[JoyMonteiro]

Just to make more explicit, you can change the arguments to the boundary layer like this

simple_physics = SimplePhysics(top_of_boundary_layer=50000.0)

[MengZuo]

Hi Joy, thank you very much for your detailed reply! And sorry for the late reply.
After modifying the surface pressure followed by your instructions, the model can run correctly, but the results are different from what I anticipated. I want to remove the physical processes below 700hPa, preventing them from being involved in the calculations, in order to simulate a RCE state in higher terrain region. Could you please advise if there are any other parameters that need to be adjusted? Thank you very much, and looking forward to your reply.

[MengZuo]

Hi Joy, thank you very much for your detailed reply! And sorry for the late reply.
After modifying the surface pressure followed by your instructions, the model can run correctly, but the results are different from what I anticipated. I want to remove the physical processes below 700hPa, preventing them from being involved in the calculations, in order to simulate a RCE state in higher terrain region. Could you please advise if there are any other parameters that need to be adjusted? Thank you very much, and looking forward to your reply.

[JoyMonteiro]

Hello Meng Zuo, sorry for the delayed reply.

Could you please send me the script that you have written, and explain a little more what you anticipated and what results you have obtained? @Ai33L might have written some components which may also be of help.

[MengZuo]

Hi Joy, thank you very much for your reply. @JoyMonteiro @Ai33L
I am trying to modify the surface air pressure. My purpose is to set the surface air pressure to 700hPa or 500hPa and analyze the elevated heating effect, for example, over the high terrain regions such as the Tibetan Plateau in Asia and the Colorado Plateau in North [America.]
The script is based on the example you given in the climt tools: gmd_radiative_convective.py
I just modified the code as is shown below.
`from sympl import (
PlotFunctionMonitor, AdamsBashforth, NetCDFMonitor
)
from climt import SimplePhysics, get_default_state, get_grid
import numpy as np
from datetime import timedelta

from climt import EmanuelConvection, RRTMGShortwave, RRTMGLongwave, SlabSurface
import matplotlib.pyplot as plt

def plot_function(fig, state):
ax = fig.add_subplot(2, 2, 1)
ax.plot(
state['air_temperature_tendency_from_convection'].to_units(
'degK day^-1').values.flatten(),
state['air_pressure'].to_units('mbar').values.flatten(), '-o')
ax.set_title('Conv. heating rate')
ax.set_xlabel('K/day')
ax.set_ylabel('millibar')
ax.grid()

ax.axes.invert_yaxis()
ax = fig.add_subplot(2, 2, 2)
ax.plot(
    state['air_temperature'].values.flatten(),
    state['air_pressure'].to_units('mbar').values.flatten(), '-o')
ax.set_title('Air temperature')
ax.axes.invert_yaxis()
ax.set_xlabel('K')
ax.grid()

ax = fig.add_subplot(2, 2, 3)
ax.plot(
    state['air_temperature_tendency_from_longwave'].values.flatten(),
    state['air_pressure'].to_units('mbar').values.flatten(), '-o',
    label='LW')
ax.plot(
    state['air_temperature_tendency_from_shortwave'].values.flatten(),
    state['air_pressure'].to_units('mbar').values.flatten(), '-o',
    label='SW')
ax.set_title('LW and SW Heating rates')
ax.legend()
ax.axes.invert_yaxis()
ax.set_xlabel('K/day')
ax.grid()
ax.set_ylabel('millibar')

ax = fig.add_subplot(2, 2, 4)
net_flux = (state['upwelling_longwave_flux_in_air'] +
            state['upwelling_shortwave_flux_in_air'] -
            state['downwelling_longwave_flux_in_air'] -
            state['downwelling_shortwave_flux_in_air'])
ax.plot(
    net_flux.values.flatten(),
    state['air_pressure_on_interface_levels'].to_units(
        'mbar').values.flatten(), '-o')
ax.set_title('Net Flux')
ax.axes.invert_yaxis()
ax.set_xlabel('W/m^2')
ax.grid()
plt.tight_layout()
plt.savefig("test_rce_ps_to_700hpa.pdf", dpi=300, format="pdf")

monitor = PlotFunctionMonitor(plot_function)

timestep = timedelta(minutes=5)

convection = EmanuelConvection()
radiation_sw = RRTMGShortwave()
radiation_lw = RRTMGLongwave()
slab = SlabSurface()
simple_physics = SimplePhysics()

simple_physics = SimplePhysics(top_of_boundary_layer=50000.0)

store_quantities = ['air_temperature',
'air_pressure',
'specific_humidity',
'air_pressure_on_interface_levels',
'air_temperature_tendency_from_convection',
'air_temperature_tendency_from_longwave',
'air_temperature_tendency_from_shortwave']
netcdf_monitor = NetCDFMonitor('rad_conv_eq_ps_to_700hpa.nc',
store_names=store_quantities,
write_on_store=True)
convection.current_time_step = timestep

state = get_default_state([simple_physics, convection,
radiation_lw, radiation_sw, slab])

grid_state = get_grid(p_surf_in_Pa=70000.0)

state = get_default_state([simple_physics, convection,
radiation_lw, radiation_sw, slab], grid_state=grid_state)

state['air_temperature'].values[:] = 270
state['surface_albedo_for_direct_shortwave'].values[:] = 0.5
state['surface_albedo_for_direct_near_infrared'].values[:] = 0.5
state['surface_albedo_for_diffuse_shortwave'].values[:] = 0.5

state['zenith_angle'].values[:] = np.pi/2.5
state['surface_temperature'].values[:] = 300.
state['ocean_mixed_layer_thickness'].values[:] = 5
state['area_type'].values[:] = 'sea'

time_stepper = AdamsBashforth([convection, radiation_lw, radiation_sw, slab])

for i in range(20000):
convection.current_time_step = timestep
diagnostics, state = time_stepper(state, timestep)
state.update(diagnostics)
diagnostics, new_state = simple_physics(state, timestep)
state.update(diagnostics)
if (i+1) % 20 == 0:
monitor.store(state)
netcdf_monitor.store(state)
print(i, state['surface_temperature'].values)
print(state['surface_upward_sensible_heat_flux'])
print(state['surface_upward_latent_heat_flux'])

state.update(new_state)
state['time'] += timestep
state['eastward_wind'].values[:] = 3.

`

And what I anticipate is to get the following result, as Hu and Boos done in their paper based on the CliMT tool.

Hu, S., and W. R. Boos, 2017: The Physics of Orographic Elevated Heating in Radiative–Convective Equilibrium. J. Atmos. Sci., 74, 2949–2965, https://doi.org/10.1175/JAS-D-16-0312.1.

Thank you very much, and looking forward to your reply.

[Ai33L]

Hi @MengZuo

Sorry for the delayed response, just seeing this now.

I ran your script for a few values of surface pressures. This is what I get, seems similar to the plot you sent -

As Joy said, setting surface pressure as an argument and setting top of boundary layer above that level should be sufficient in your case. There should also be no processes below set surface pressure that come into calculations.

I find that the model needs around 6 months to equilibriate from initialisation. The results can also vary a bit with surface albedo and how you set boundary layer height.

Hope this helps! Please let me know if there's something I missed out.

[MengZuo]

Hi @Ai33L @JoyMonteiro
Thank you very much for your reply.
When I set the runtime to 6 months, the results seemed reasonable. Thank you!
By the way, there is another question. The default area type in the model is set to 'sea'. To represent the reduced
evaporation that typically occurs over land compared to ocean, I want to multiply the latent heat flux by a constant
fraction 0.5 or 0.25.
I noticed that the area type can be set to 'land', state['area_type'].values[:] = 'land', but I'm not clear on which specific parameters change after adjusting sea to land? Looking forward to your reply, thank you!

[Ai33L]

Hi @MengZuo

In your model configuration, area type sets suitable default values for heat capacity, depth and surface material density in the SlabSurface component. You can find the relevant code here -

climt/climt/_components/slab_surface.py

Lines 113 to 125 in c3bdf05

	land_mask = np.logical_or(area_type == 'land', area_type == 'land_ice')
	sea_mask = np.logical_or(area_type == 'sea', area_type == 'sea_ice')
	land_ice_mask = area_type == 'land_ice'
	sea_ice_mask = area_type == 'sea_ice'
	net_heat_flux[land_ice_mask] = -raw_state['upward_heat_flux_at_ground_level_in_soil'][land_ice_mask]
	net_heat_flux[sea_ice_mask] = raw_state['heat_flux_into_sea_water_due_to_sea_ice'][sea_ice_mask]
	raw_state['surface_material_density'][sea_mask] = raw_state['sea_water_density'][sea_mask]
	raw_state['surface_thermal_capacity'][land_mask] = raw_state['heat_capacity_of_soil'][land_mask]
	diagnostics['depth_of_slab_surface'][sea_mask] =\
	raw_state['ocean_mixed_layer_thickness'][sea_mask]
	diagnostics['depth_of_slab_surface'][land_mask] =\
	raw_state['soil_layer_thickness'][land_mask]

The easiest way to implement reduced latent heat flux over land in the Simple physics component is to manually set a lower value for surface specific humidity. You can do this by making these changes to your code -

add

from climt import bolton_q_sat
from sympl import get_constant

change

simple_physics = SimplePhysics(use_external_surface_specific_humidity=True, top_of_boundary_layer=50000.0)

use_external_surface_specific_humidity allows you to manually set surface specific humidity. You've declared simple_physics twice in your code, just do it once.

and add these four lines to calculate surface specific humidity before you call simple_physics (inside the time loop). You can change scaling as you wish.

scaling=0.5
Rd = get_constant('gas_constant_of_dry_air', 'J kg^-1 K^-1')
Rh2o = get_constant('gas_constant_of_vapor_phase', 'J/kg/degK')
state['surface_specific_humidity'][:]= bolton_q_sat(state['surface_temperature'], state['surface_air_pressure'], Rd, Rh2o) * scaling
    
diagnostics, new_state = simple_physics(state, timestep)
state.update(diagnostics)

Hope this helps!