Goal: In this example, I want to give the Meteorology component a constant scalar precipitation value and check whether it produces output when the model state is updated.
Start with an import and some magic:
In [ ]:
%matplotlib inline
import numpy as np
Import the Meteorology component and create an instance:
In [ ]:
from topoflow.components.met_base import met_component
m = met_component()
Locate the cfg file and initialize the component:
In [ ]:
cfg_file = './input/meteorology.cfg'
m.initialize(cfg_file)
Despite setting a value of 20 mm/hr for P in the cfg file, if I call get_value at this point, the model precip values are zeroed out:
In [ ]:
precip = m.get_value('atmosphere_water__precipitation_leq-volume_flux') # `P` internally
print type(precip)
print precip.size
precip
Maybe this will change after the first update?
Show the model start, current, and stop times:
In [ ]:
print m.get_start_time(), m.get_current_time(), m.get_end_time()
Advance the model by one time step:
In [ ]:
m.update()
print '\nCurrent time: {} s'.format(m.get_current_time())
Note that it hasn't precipitated:
In [ ]:
precip = m.get_value('atmosphere_water__precipitation_leq-volume_flux')
print precip
But this might be expected, since the initial precip value retrieved by get_value was also zero.
Try manually setting the precip rate variable to an arbitrary non-zero value:
In [ ]:
new_value = np.array(15, dtype=np.float64) # set_value doesn't convert to the TF type...
m.set_value('atmosphere_water__precipitation_leq-volume_flux', new_value)
precip = m.get_value('atmosphere_water__precipitation_leq-volume_flux')
print precip
Note that I can't simply give set_value a new value of 15; I need to create a one-element Numpy array with a value of 15 to match what TopoFlow needs internally.
Advance the model by another time step:
In [ ]:
m.update()
print '\nCurrent time: {} s'.format(m.get_current_time())
Is precipitation being produced?
In [ ]:
precip = m.get_value('atmosphere_water__precipitation_leq-volume_flux')
print precip
No, it's not.
Look at the verbose output from the call to update:
Calling update_P_integral()...
Calling update_P_max()...
Calling update_P_rain()...
>> Rain is falling...
Calling update_P_snow()...
Calling update_bulk_richardson_number()...
Calling update_bulk_aero_conductance()...
Callilng update_sensible_heat_flux()...
Calling update_saturation_vapor_pressure()...
Calling update_saturation_vapor_pressure()...
Calling update_vapor_pressure()...
Calling update_dew_point()...
Calling update_precipitable_water_content()...
Calling update_vapor_pressure()...
Calling update_latent_heat_flux()...
Calling update_conduction_heat_flux()...
Calling update_advection_heat_flux()...
Calling update_julian_day()...
Calling update_net_shortwave_radiation()...
Calling update_em_air()...
Calling update_net_longwave_radiation()...
Calling update_net_energy_flux()...
Calling read_input_files()...
Reached end of scalar rainfall duration.
P set to 0 by read_input_files().
Calling write_output_files()...The fourth line says it's raining, but the second-to-last line looks problematic.
On examining the TopoFlow source code, the following code block starts at line 1827 in met_base.py:
if (self.P_type.lower() != 'scalar'):
#------------------------------------
# Precip is unique in this respect.
#--------------------------------------------------
# 2/7/13. Note that we don't change P from grid
# to scalar since that could cause trouble for
# other comps that use P, so we just zero it out.
#--------------------------------------------------
self.P.fill( 0 )
if (self.DEBUG):
print 'Reached end of file:', self.P_file
print ' P set to 0 by read_input_files().'
So if the precipitation rate is a scalar, it's explicitly zeroed out, by design, after every update.
Result: Rain is falling, but it's canceled at the end of every update.