Additional material: Python code of the coupling prototype
The coupling code needs to have CMF and LandscapeDNDC installed as python libraries in order to run the coupled model.
# -*- coding=utf-8 -*-
from __future__ import division
import cmf as cmf
import numpy as np
import sys
from time import time
import os
class Cell(object):
'''A "bridge" object between a cmf cell and a LandscapeDNDC instance (grid)
This class combines a cmf cell with a LandscapeDNDC cell from LandscapeDNDC
The layers of the cmf cell will be created using the data of the LandscapeDNDC cell'''
def __init__(self,cmf_cell,mobile):
'''
cmf_cell : A cell of cmf (with area and location)
mobile : LandscapeDNDC instance for the same location
'''
# ensure cmf cell has no layers
if len(cmf_cell.layers):
raise RuntimeError('''Creating a combined / cmf cell requires cells without layers.
The layers will be created using the data from LandscapeDNDC''')
self.cmf = cmf_cell
self.mobile = mobile
# just a shortcut to site
site = self.mobile.site
# Gets a Brooks-Corey b (curve shape parameter) from clay content. Function was derived by fitting
# a boltzman distribution function into the Brooks-Corey Parameter estimation routine of cmf
# using Bodenkundliche Kartieranleitung 4
b = 12.958/(1+np.exp(-(np.asarray(site.clay_sl)-0.32)/0.1124)) + 1.99
depth=0.0# current soildepth of cmf cell
# For each thickness of a layer
for i,h in enumerate(self.mobile.site.h_sl):
if h>0 : # If layer exists
# Increase soil depth
depth += h
# Get porosity
poro = site.poro_sl[i]
# Get field capacity
fc = site.wcMax_sl[i] * 1e-3# convert from l/m3 to m3/m3
# Get wilting point (not used)
wilt = site.wcMin_sl[i] * 1e-3# convert from l/m3 to m3/m3
# Get saturated conductivity
Ks = site.sks_sl[i] * 24 * 60 / 100# convert from cm/min to m/day
# Create a Brooks-Corey retention curve for the data
bc = cmf.BrooksCoreyRetentionCurve(Ks,poro,b[i],fc)
# Create the layer
self.cmf.add_layer(depth, bc)
# Set the layer's water content
self.cmf.layers[-1].theta = self.mobile.watercycle.wc_sl[i] * 1e-3# convert from l/m3 to m3/m3
print"%s: Added %i layers, soildepth is %g m" % (self.cmf,len(self.cmf.layers),self.cmf.soildepth)
self.nlayers = self.cmf.layer_count()
self.cumN2O = 0.0
def __repr__(self):
return"<LandscapeDNDC/cmf Cell %s,%s>" % (self.mobile, self.cmf)
def update_cmf(self):
"""Gets data from LandscapeDNDC to update parameters of cmf
Parameters updated by LandscapeDNDC only (LandscapeDNDC->cmf):
LAI
Vegetation height
concentration of throughfall
Parameters updated by LandscapeDNDC, to be changed by cmf (LandscapeDNDC<->cmf)
Solute state of each layer and each solute
"""
# Update vegetation (self.cmf.Vegetation.LAI = xxx, fails without notification, bug to be handled in cmf)
v=self.cmf.vegetation
v.LAI = max(self.mobile.physiology.lai_vt.sum(),1)
v.Height = max(self.mobile.site.hMax_vt)
v.StomatalResistance=200
self.cmf.vegetation=v
# Get soilchemistry
sc = self.mobile.soilchemistry
ph = self.mobile.physiology
mc = self.mobile.microclimate
A = self.cmf.area # Area of the cell
# Convert kg m-2 to g (for states and fluxes)
sc_states = (sc.nh4_sl,sc.nh3_sl,sc.no3_sl,sc.urea_sl,sc.don_sl,sc.doc_sl)
A = self.cmf.area # Area of the cell
# Set solute storages and fluxes
for j,X in enumerate(self.cmf.project.solutes):
for i,l in enumerate(self.cmf.layers):
ifnot np.isfinite(sc_states[j][i+1]) and np.isfinite(l.Solute(X).state):
print > sys.stderr, "LandscapeDNDC created infinite value for %s in %s \n CMF is not updated and will overwrite the value in LandscapeDNDC" % (X,l)
print > sys.stderr, sc_states[j][i+1]
else:
# Convert from kg/m2 -> g
l.Solute(X).state = sc_states[j][i+1] * A * 1e3
# Put throughfall into first layer if the ground is dry, else into surfacewater
throughfall_target = self.cmf.layers[0] #if self.cmf.surfacewater.is_empty() else self.surfacewater
# Convert kg m-2 to g (for states and fluxes) = 1e3 * A
sc_fluxes = 1e3 * A * np.array((ph.throu_nh4,0.0,ph.throu_no3))
# Set throughfall flux source in first layer or surface water
# TODO: Was passiert eigentlich mit den fluxes?
for j,X in enumerate(list(self.cmf.project.solutes)[:3]):
throughfall_target.Solute(X).source = sc_fluxes[j]
# Get current weather conditions from LandscapeDNDC
w=cmf.Weather(mc.temp_nd,mc.tMin_nd,mc.tMax_nd,cmf.rH_from_vpd(mc.temp_nd,mc.vpd_nd),mc.win_nd,0.5,mc.rad_nd * (24*60*60*1e-6))
w.instrument_height=2
self.cmf.set_weather(w)
self.cmf.set_rainfall(mc.prec_nd)
# integrate fluxes
self.cumN2O += self.mobile.soilchemistry.nd_dN2O
@property
def biomass(self):
ph=self.mobile.physiology
return ph.mFol_vt.sum() + ph.mFrt_vt.sum() + ph.mBud_vt.sum() + ph.mSap_vt.sum()
def update_mobile(self,t,coupled):
"""Gets data from cmf to update parameters of LandscapeDNDC
Parameters updated by cmf only (cmf->LandscapeDNDC):
Water content for each layer
Parameters updated by cmf, to be changed by LandscapeDNDC (LandscapeDNDC<->cmf)
Solute state of each layer and each solute
"""
# Update water content in mobile.watercycle
m=self.mobile
wc = m.watercycle
#sc= m.soilchemistry
layers = self.cmf.layers[0:self.nlayers]
wc.wc_sl[1:self.nlayers+1] = layers.wetness * layers.porosity * 1e3
m.site.WTDEPTH = self.cmf.saturated_depth
wc.transp = self.cmf.transpiration.waterbalance(t)
wc.evacep = self.cmf.canopy.flux_to(self.cmf.evaporation,t) if self.cmf.canopy else0.0
wc.evasoil = self.cmf.layers[0].flux_to(self.cmf.evaporation,t)
wc.evapot = cmf.PenmanMonteith(self.cmf.get_weather(t),self.cmf.vegetation,2.0)
if coupled :
# Get soilchemistry
sc=self.mobile.soilchemistry
sc_states = (sc.nh4_sl,sc.nh3_sl,sc.no3_sl,sc.urea_sl,sc.don_sl,sc.doc_sl)
A = self.cmf.area # Area of the cell
# Set solute storages
for j,X in enumerate(self.cmf.project.solutes):
for i,l in enumerate(self.cmf.layers):
if np.isfinite(sc_states[j][i+1]) andnot np.isfinite(l.Solute(X).state):
print > sys.stderr, "cmf created infinite value for %s in %s LandscapeDNDC is not updated and will overwrite the value in cmf" % (X,l)
else:
sc_states[j][i+1] = l.Solute(X).state * 1e-3 / A
#endif coupled
def make_outlet(cell):
# Outlet on left hand side of lowest cell
opotential = 2.9
ocell = 13
outlet = p.NewOutlet("Outlet",Cells[ocell].cmf.x-10,Cells[ocell].cmf.y,Cells[ocell].cmf.z-opotential)
outlet.potential = Cells[ocell].cmf.z - opotential
# Connect the outlet with the left hand cell
return outlet
def make_solver(project):
# Create a solver for the water volume ODE's
w_solver = cmf.CVodeIntegrator(1e-6)
w_solver.preconditioner = 'R'
# Create a solver for solute transport ODE's
s_solver = cmf.ImplicitEuler(1e-6)
# Combine solvers
solver = cmf.SoluteWaterIntegrator(w_solver,s_solver,p)
solver(cmf.h)
return solver
def run(cmf_solver, M2d, Cells,profiling=False,coupled=True):
"""Should demonstrate how the coupled system will be solved for ONE timestep of length dt
"""
# get timestep length from m2d
m2d_timestep = cmf.sec * ((24 * 3600) / M2d.time.tsMax())
tstart = time()
M2d.launch()
# set LandscapeDNDC init NO3 concentrations to zero
tM2d = time()
for C in Cells:
C.update_cmf()
t_update_cmf=time()
end = cmf_solver.t + m2d_timestep
while cmf_solver.t<end:
cmf_solver(cmf.h)
t_cmf = time()
for C in Cells:
C.update_mobile(cmf_solver.t,coupled)
t_update_Mobile=time()
if profiling:
print"%-24s (Total:%6.2fs, M2d:%5.2fs, cmf:%5.2fs)" % (cmf_solver.t,t_update_Mobile-tstart,tM2d-tstart,t_cmf-t_update_cmf)
if"plot"in sys.argv or"final"in sys.argv:
import matplotlib.pyplot as plt
class plotter(object):
"""A plotter for a cmf/LandscapeDNDC hillslope
Creates two subplots, one for water content and one for a solute (NO3 by default)
"""
def __call__(self):
self.wet(solver.t)
for c,bb,bn in zip(Cells,self.biomassbars,self.N2Obars):
bb.set_height(c.biomass * 3)
bn.set_height(c.cumN2O * 2e4)
self.solute(solver.t)
def __init__(self,solute):
plt.subplot(211)
self.wet = cmf.draw.hill_plot(p,solver.t,cmap=plt.cm.jet_r)
self.wet.scale = 20
plt.xlabel(" ")
plt.ylabel("height [m]")
plt.ylim(-3,14)
plt.xlim(0,260)
plt.subplot(212)
self.solute = cmf.draw.hill_plot(p,solver.t,solute,cmap=plt.cm.jet)
self.solute.scale = 1e100
self.solute.evalfunction = lambda l: (l.conc(solute) * 0.01) **0.25
plt.ylabel("height [m]")
plt.xlabel("distance [m]")
font = {'size' : '6.0'}
plt.rc('font', **font)
plt.rc('lines', linewidth=0.5, color='b')
left=np.array([c.cmf.x for c in Cells])
bottom=np.array([c.cmf.z for c in Cells])
self.biomassbars = plt.bar(left+0-2,1e-2*np.ones(len(Cells)),
width=4.0,bottom=bottom,color='g',zorder=20, linewidth=0)
self.N2Obars = plt.bar(left+0+2,1e-2*np.ones(len(Cells)),
width=4.0,bottom=bottom,color='r',zorder=20, linewidth=0)
plt.ylim(-3,14)
plt.xlim(0,260)
self()
def set_scale(self,new_scale):
self.solute.scale=new_scale
self.wet.scale=new_scale
self()
if __name__=='__main__':
print cmf.VERSION
# Interesting code begins...
print"loading LandscapeDNDC"
import m2d
print"loading cmf"
import cmf
# Create a log file
m2d.M2dLog.init('log.txt',5)
# Create the LandscapeDNDC object using the input files
M2d = m2d.Mobile2D(m2d.OPT_NONE,"input2D/setup2D.xml")
print" -> LandscapeDNDC object generated."
print"Create a new CMF project with solutes: NH4 NH3 NO3 urea DON DOC "
p = cmf.project('NH4 NH3 NO3 urea DON DOC')
NH4, NH3, NO3, urea, DON, DOC = p.solutes
Cells = []
slope = 0.5
height = 0.1
offset = 0
area = 100.0
for i,m in enumerate(M2d):
# Create a hillslope with cmf cell with 100m2 area
if i ==0:
z = (12-i) * slope * 0.5 + height + 5 * slope * 0.5
y = 0.0
x = 10.0 * (i-offset)
c=p.NewCell(x, y, z, area)
elif i < 8 :
z = (12-i) * slope * 0.5 + height + 5 * slope * 0.5
y = 0.0
x = 10.0 * (i-offset)
c=p.NewCell(x, y, z, area)
elif (i < 13) and (i>=8) :
y = 0.0
z = (12-i) * slope + height
x = 10.0 * (i-offset)
c=p.NewCell(x, y, z, area)
elif i == 12 :
y = 0.0
z = height
x = 10 * (i-offset)
c=p.NewCell(x,y, z, area)
elif i == 13 :
z = 0
y = 0.0
x= 10 * (i-offset)
c=p.NewCell(x, y, z, area)
elif i == 14 :
z = height
x = 10 * (i-offset)
y = 0.0
c=p.NewCell(x,y , z, area)
elif (i >= 15) and (i < 19) :
x = 10.0 * (i-offset)
y = 0.0
z = height + slope * (i-14)
c=p.NewCell(x, y, z, area)
elif (i >=19) and (i<26) :
x = 10.0 * (i-offset)
y = 0.0
z = height + 0.5 * slope * (i-14) + 5 * slope * 0.5
c=p.NewCell(x, y, z, area)
else :
x = 10.0 * (i-offset)
y = 0.0
z = height + 0.5 * slope * (i-14) + 5 * slope * 0.5
c=p.NewCell(x, y, z, area )
# Mesh the hillslope (catena, 1d chain topology)
if i: c.topology.AddNeighbor(p[-2],10)
# Create a combination cell
Cells.append(Cell(c,m))
# Use Richard's equation for percolation
c.install_connection(cmf.Richards)
# Use Penman-Monteith eq. ET
#endfor
print" -> CMF created 2d hill slope. "
# cmf setup stuff, making lateral flow and an outlet at foot of hillslope
cmf.connect_cells_with_flux(p.cells,cmf.Richards_lateral)
outlet = make_outlet(p[13])
solver=make_solver(p)
# Get the start date from LandscapeDNDC
DOY, year = M2d.time.nd()-1, M2d.time.ny()
solver.t = cmf.Time(1,1,year) + DOY * cmf.day
C=Cells[0] # Shortcut for tests
print"Let us start..."
# If "plot" is in the arguments, create a class for simple plotting
if"plot"in sys.argv:
if"show"in sys.argv:
plt.figure(num=None, figsize=(4, 4.5), dpi=300, facecolor='w', edgecolor='k')
hp = plotter(NO3)
plt.show()
plt.clf()
#endif
plt.figure(num=None, figsize=(4, 4.5), dpi=300, facecolor='w', edgecolor='k')
hp = plotter(NO3)
plt.ylabel("height [m]")
plt.ioff()
# Create picture directory if it is not there
ifnot os.path.exists('pic'):
os.mkdir('pic')
#endif
# get all layers
all_layers=cmf.node_list(np.ravel([c.layers for c in p]))
#endif
days = 10*365 + 1
coupled = True
if"uncoupled"in sys.argv:
coupled = False
if"run"in sys.argv:
for i in range(1,days):
run(solver,M2d,Cells,True,coupled)
if"plot"in sys.argv and i in range (1,days,1):
hp()
plt.draw()
plt.savefig('pic/coupled_%04i.pdf' % i, dpi=300, format="pdf", transparent=True)
print"Max cum dN2O=%g kg/ha" % max(c.cumN2O for c in Cells)
#endif
# reset the yearly accumulation
if365==i:
for c in Cells :
c.cumN2O = 0
if730==i:
for c in Cells :
c.cumN2O = 0
if1095==i:
for c in Cells :
c.cumN2O = 0
if1461==i:
for c in Cells :
c.cumN2O = 0
#endfor
#endif
if"final"in sys.argv :
hp = plotter(NO3)
hp()
plt.show()
#endif
#end