CliMT / climt

The official home of climt, a Python based climate modelling toolkit.
https://climt.readthedocs.io
Other
160 stars 46 forks source link

RRTMG SW behaviour undocumented #112

Open suhasdl opened 5 years ago

suhasdl commented 5 years ago

Description

I was trying to vary the solar constant by changing the stellar_irradiance value. However when its set to 0 the solar constant defaults to the original value of about 1360. But it works fine even for a small value of say 10^-3.

I am attaching a sample code which demonstrates this.

#########################################
import climt
from sympl import (
    PlotFunctionMonitor, NetCDFMonitor,
    TimeDifferencingWrapper, UpdateFrequencyWrapper,
    set_constant, DataArray
)
import numpy as np
from datetime import timedelta

set_constant('stellar_irradiance', value=0, units='W m^-2')
model_time_step = timedelta(minutes=10)

# Create components
convection     = climt.EmanuelConvection()
simple_physics = TimeDifferencingWrapper(climt.SimplePhysics())
radiation_step = timedelta(hours=1)
radiation_lw   = UpdateFrequencyWrapper(
    climt.RRTMGLongwave(), radiation_step)
radiation_sw   = UpdateFrequencyWrapper(
    climt.RRTMGShortwave(), radiation_step)
slab_surface   = climt.SlabSurface()
dycore = climt.GFSDynamicalCore(
    [simple_physics, slab_surface, radiation_sw,
     radiation_lw, convection], number_of_damped_levels=5)

grid = climt.get_grid(nx=128, ny=64)
# Create model state
my_state = climt.get_default_state([dycore], grid_state=grid)

# Set initial/boundary conditions
latitudes  = my_state['latitude'].values
longitudes = my_state['longitude'].values
surface_shape = latitudes.shape
zenith_angle = np.radians(latitudes)
my_state['zenith_angle'].values = zenith_angle

# Set SST profile
my_state['surface_temperature'] = DataArray(
    300.*np.ones(surface_shape),
    dims=['lat', 'lon'], attrs={'units': 'degK'})

my_state['eastward_wind'].values[:] = np.random.randn(
    *my_state['eastward_wind'].shape)

my_state['ocean_mixed_layer_thickness'].values[:] = 50

tmax  = my_state['time'] + timedelta(days=1000)
tNext = my_state['time']

while (my_state['time'] <= tmax):
    diag, my_state = dycore(my_state, model_time_step)
    my_state.update(diag)
    my_state['time'] += model_time_step
    print('max. downwelling short wave flux: ', np.amax(my_state['downwelling_shortwave_flux_in_air'].values))

#####################################################
JoyMonteiro commented 5 years ago

This is simply the behaviour of the RRTMG SW code. See the comments in https://github.com/CliMT/climt/blob/69d1e5d7071f248da5c68053ce08af2a11267cd1/climt/_lib/rrtmg_sw/rrtmg_sw_rad.f90#L269-L279

This should be documented better. I will keep this issue open to remind me of the same.

suhasdl commented 5 years ago

Thanks Joy.