Precipitation in the Meteorology component

Goal: In this example, I give the Meteorology component a time series of precipitation values and check whether it produces output when the model state is updated.

Define a helpful constant:



In [ ]:

    
mps_to_mmph = 1000 * 3600

Programmatically create a file holding the precipitation rate time series. This will mimic what I'll need to do in WMT, where I'll have access to the model time step and run duration. Start by defining the precipitation rate values:



In [ ]:

    
import numpy as np
n_steps = 10  # can get from cfg file
precip_rates = np.linspace(5, 20, num=n_steps, endpoint=False)
precip_rates

Next, write the values to a file to the input directory, where it's expected by the cfg file:



In [ ]:

    
np.savetxt('./input/precip_rates.txt', precip_rates, fmt='%6.2f')

Check the file:



In [ ]:

    
cat input/precip_rates.txt

BMI component

Import the BMI Meteorology component and create an instance:



In [ ]:

    
from topoflow.components.met_base import met_component
m = met_component()

Initialize the model. A value of snow depth h_snow is needed for the model to update.



In [ ]:

    
m.initialize('./input/meteorology-1.cfg')
m.h_snow = 0.0  # Needed for update

Unlike when P is a scalar, the initial model precipitation volume flux is the first value from precip_rates.txt:



In [ ]:

    
precip = m.get_value('atmosphere_water__precipitation_leq-volume_flux')  # `P` internally
print type(precip)
print precip.size
precip * mps_to_mmph

Advance the model by one time step:



In [ ]:

    
m.update()
print '\nCurrent time: {} s'.format(m.get_current_time())

Unlike the scalar case, there's an output volume flux of precipitation:



In [ ]:

    
print precip * mps_to_mmph  # note that this is a reference, so it'll take the current value of `P`

Advance the model to the end, saving the model time and output P values (converted back to mm/hr for convenience) at each step:



In [ ]:

    
time = [m.get_current_time().copy()]
flux = [precip.copy() * mps_to_mmph]
while m.get_current_time() < m.get_end_time():
    m.update()
    time.append(m.get_current_time().copy())
    flux.append(m.get_value('atmosphere_water__precipitation_leq-volume_flux').copy() * mps_to_mmph)

Check the time and flux values:



In [ ]:

    
time



In [ ]:

    
flux

Result: Fails. Input precipipation rates do not match output precipitation volume flux because of changes we made to TopoFlow source.

Babel-wrapped component

Import the Babel-wrapped Meteorology component and create an instance:



In [ ]:

    
from cmt.components import Meteorology
met = Meteorology()

Initialize the model.



In [ ]:

    
%cd input
met.initialize('meteorology-1.cfg')

The initial model precipitation volume flux is the first value from precip_rates.txt:



In [ ]:

    
bprecip = met.get_value('atmosphere_water__precipitation_leq-volume_flux')
print type(bprecip)
print bprecip.size
print bprecip.shape
bprecip * mps_to_mmph

Advance the model to the end, saving the model time and output P values (converted back to mm/hr for convenience) at each step:



In [ ]:

    
time = [met.get_current_time()]
flux = [bprecip.max() * mps_to_mmph]
count = 1
while met.get_current_time() < met.get_end_time():
    met.update(met.get_time_step()*count)
    time.append(met.get_current_time())
    flux.append(met.get_value('atmosphere_water__precipitation_leq-volume_flux').max() * mps_to_mmph)
    count += 1

Check the time and flux values (noting that I've included the time = 0.0 value here):



In [ ]:

    
time



In [ ]:

    
flux

Result: Fails. Input precipipation rates do not match output precipitation volume flux because of changes we made to TopoFlow source.