How to specify initial condition for column RCE model

[spencerahill]

Searching through the docs and source code, I can't figure out how to specify initial conditions for an RCE model (either a single column or lat-by-lat). Motivation: I am exploring cases with very high albedos, leading to very cold equilibrium climates. But then it takes a long time (many hundreds of days) to reach steady state when starting at the near-Earth-like temperatures of the default initial condition.

Can the initial condition be specified, and if so how? TIA.

[brian-rose]

This one is easy: just slice the state variable arrays.

For example, this creates an isothermal initial condition:

import climlab
state = climlab.column_state(num_lev=20)
state['Tatm'][:] = state['Ts']
rad = climlab.radiation.RRTMG(state=state)
print(rad.state)

[brian-rose]

This relates to #49

The docs could use a major overhaul to be more user-centric. I envision a "climlab recipes" section full of simple and not-so-simple working examples.

[brian-rose]

You can accomplish the same thing by assembling a model and then slicing any of the state variable arrays, since they all point to the same memory.

Just for fun, here is an example of an RCM with interactive water vapor and surface turbulent fluxes. Start it from a warm and dry initial condition and watch it moisten itself and adjust to a convectively neutral lapse rate!

import climlab
from climlab import constants as const
short_timestep = const.seconds_per_hour
long_timestep = short_timestep*6
full_state = climlab.column_state(water_depth=10.)
#  Initialize a nearly dry column (small background stratospheric humidity)
full_state['q'] = 0.*full_state.Tatm + 5.E-6
temperature_state = {'Tatm':full_state.Tatm,'Ts':full_state.Ts}
#  Surface model
shf = climlab.surface.SensibleHeatFlux(name='Sensible Heat Flux',
                           state=temperature_state, Cd=0.5E-3,
                           timestep=short_timestep)
lhf = climlab.surface.LatentHeatFlux(name='Latent Heat Flux',
                         state=full_state, Cd=0.5E-3,
                         timestep=short_timestep)
surface = climlab.couple([shf,lhf], name="Slab")
#  Convection scheme -- water vapor is a state variable
conv = climlab.convection.EmanuelConvection(name='Convection',
                             state=full_state,
                             IPBL=1,
                             timestep=short_timestep,)
rad = climlab.radiation.RRTMG(name='Radiation',
                    state=temperature_state,
                    specific_humidity=full_state.q,
                    albedo=0.18,
                    insolation=341.3,
                    timestep=long_timestep,
                    )
atm = climlab.couple([rad, conv], name='Atmosphere')
model = climlab.couple([atm,surface], name='Moist column model')
#  Set a warm isothermal initial condition
initial_temp = 300.
model.state['Ts'][:] = initial_temp
model.state['Tatm'][:] = initial_temp
print(model)

[spencerahill]

Thanks Brian, that's exactly what I needed. I agree re: clarifying this in the docs, as the solution is ultimately quite simple. (You might even consider just using one or both of the examples you gave above for that purpose.)

Separately, it might ultimately be nice to support specifying the ICs via a kwarg in the state constructor. Even more user-friendly than having to operate on items in the state object's dictionary.

I consider this closed, but your call given the docs considerations.

[spencerahill]

Another recipe re: initial conditions worth considering for the docs is to initialize it with the radiative equilibrium solution from a simple 1 layer greenhouse model, on the assumption that the RCE solution will probably be not too far from whatever the local rad-eq solution is, given the insolation and albedo. E.g. something like

STEF_BOLTZ_CONST = 5.6704e-8
TEMP_MIN_VALID = 150.
TEMP_MAX_VALID = 340.

# Leading 2^(1/4) power is from simple 1 layer greenhouse.
temp_rad_eq = (2**0.25)*(insol*(1-albedo) / STEF_BOLTZ_CONST)**0.25

# Prevent runaway cold state or runaway greenhouse (this part is probably optional re: docs)
temp_rad_eq = max(temp_rad_eq, TEMP_MIN_VALID)
temp_rad_eq = min(temp_rad_eq, TEMP_MAX_VALID)

state['Ts'][:] = temp_rad_eq
state['Tatm'][:] = state['Ts'][:]

Moreover, since an isothermal atmosphere is highly stable, one could introduce a non-zero lapse rate, so that the radiation doesn't have to fight against this stable stratification for a long time before the convection can kick in.

[brian-rose]

I like these ideas. I think that a set of keyword arguments like 'isothermal', 'radiative_equilibrium', 'dry_adiabat', 'moist_adiabat' could be implemented for the column_state() function.

[brian-rose]

And another obvious one: 'standard_atmosphere'

[spencerahill]

To compute lapse rates, you need the pressure values. But as far as I can tell, the state objects don't have them:

In [6]: import climlab, pdir

In [7]: pdir(climlab.column_state())
Out[7]:
property:
    _MutableMapping__marker, _abc_impl, _allow_invalid_attributes
special attribute:
    __class__, __dict__, __doc__, __module__, __slots__, __weakref__
abstract class:
    __abstractmethods__, __subclasshook__
arithmetic:
    __add__, __radd__
object customization:
    __format__, __hash__, __init__, __new__, __repr__, __sizeof__, __str__
rich comparison:
    __eq__, __ge__, __gt__, __le__, __lt__, __ne__
attribute access:
    __delattr__, __dir__, __getattr__, __getattribute__, __setattr__
class customization:
    __init_subclass__
container:
    __contains__, __delitem__, __getitem__, __iter__, __len__, __reversed__, __setitem__
pickle:
    __getstate__, __reduce__, __reduce_ex__, __setstate__
descriptor:
    _constructor: class classmethod with getter, classmethod(function) -> method
    _valid_name: class classmethod with getter, classmethod(function) -> method
    fromkeys: class classmethod_descriptor with getter, Create a new dictionary with keys from iterable and values set to value.
class:
    _sequence_type: Built-in immutable sequence.
function:
    _build: Conditionally convert an object to allow for recursive mapping
    _configuration: The configuration for an attrmap instance.
    _delattr: Delete an attribute from the object, without attempting to
    _setattr: Add an attribute to the object, without attempting to add it as
    clear: D.clear() -> None.  Remove all items from D.
    copy: D.copy() -> a shallow copy of D
    get: Return the value for key if key is in the dictionary, else default.
    items: D.items() -> a set-like object providing a view on D's items
    keys: D.keys() -> a set-like object providing a view on D's keys
    pop: D.pop(k[,d]) -> v, remove specified key and return the corresponding value.
    popitem: D.popitem() -> (k, v), remove and return some (key, value) pair as a
    setdefault: Insert key with a value of default if key is not in the dictionary.
    update: D.update([E, ]**F) -> None.  Update D from dict/iterable E and F.
    values: D.values() -> an object providing a view on D's values
magic:
    __call__: Dynamically access a key-value pair.

In [8]: climlab.column_state()
Out[8]:
AttrDict({'Ts': Field([288.]), 'Tatm': Field([200.        , 202.68965517, 205.37931034, 208.06896552,
       210.75862069, 213.44827586, 216.13793103, 218.82758621,
       221.51724138, 224.20689655, 226.89655172, 229.5862069 ,
       232.27586207, 234.96551724, 237.65517241, 240.34482759,
       243.03448276, 245.72413793, 248.4137931 , 251.10344828,
       253.79310345, 256.48275862, 259.17241379, 261.86206897,
       264.55172414, 267.24137931, 269.93103448, 272.62068966,
       275.31034483, 278.        ])})

[brian-rose]

They actually do. The state object is just a dictionary holding objects of type Field, which is a climlab-specific hack to attach grid information to numpy arrays.

The interface to the grid information is terrible, but you can get a sense of how it works by looking at the code in climlab/domain/xarray.py that implements the translation to xarray.DataArray.

Or just do

climlab.to_xarray(state)

to access the grid information in a more user-friendly way.

[brian-rose]

I designed the Field object before I became aware of xarray data structures. The longer-term vision in #50 is to get rid of the Field object entirely.

[brian-rose]

But to be more specific: the pressure levels are in

state['Tatm'].domain.axes['lev'].points

[spencerahill]

Ah, got it. Thanks!