Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Minor updates to allometry functional unit tests #656

Merged
merged 3 commits into from
Jun 22, 2020
Merged

Minor updates to allometry functional unit tests #656

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2d57be3
Updated python allometry FUTs to py3
rgknox May 13, 2020
6bf9637
Merge branch 'master' into allometry-funcutests
rgknox May 13, 2020
f59805e
minor updates to allometry func test
rgknox May 15, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions biogeochem/FatesAllometryMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1349,6 +1349,9 @@ subroutine d2h_chave2014(d,p1,p2,p3,dbh_maxh,h,dhdd)
! Chave et al. Improved allometric models to estimate the abovegroud
! biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2015.
!
! p1 = 0.893 - E
! p2 = 0.76
! p3 = -0.034
! =========================================================================

real(r8),intent(in) :: d ! plant diameter [cm]
Expand Down Expand Up @@ -1581,6 +1584,9 @@ subroutine dh2bagw_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bagw,dbagwdd)
! Output:
! bagw: Total above ground biomass [kgC]
!
! Chave's Paper has p1 = 0.0673, p2 = 0.976
!

! =========================================================================


Expand Down
2 changes: 1 addition & 1 deletion functional_unit_testing/allometry/AutoGenVarCon.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ def __init__(self,var_sym,n_dims):


def CheckFile(filename,check_str):
file_ptr = file(filename)
file_ptr = open(filename,'r')
var_list = []
found = False
for line in file_ptr:
Expand Down
113 changes: 64 additions & 49 deletions functional_unit_testing/allometry/drive_allomtests.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,15 @@
import matplotlib.pyplot as plt
import matplotlib as mp
import ctypes
import importlib
from ctypes import * #byref, cdll, c_int, c_double, c_char_p, c_long
import xml.etree.ElementTree as ET
import argparse
import re # This is a heftier string parser
import code # For development: code.interact(local=dict(globals(), **locals()))

import sys
sys.path.append('../shared/py_src')
from PyF90Utils import c8, ci, cchar, c8_arr, ci_arr

# =======================================================================================
# Set some constants. If they are used as constant arguments to the F90 routines,
Expand Down Expand Up @@ -81,7 +84,7 @@ def setval(self,val,ipft):
def DiscreteCubeHelix(N):

base = plt.cm.get_cmap('cubehelix')
np.random.seed(1)
np.random.seed(2)
color_list = base(np.random.randint(0,high=255,size=N))
cmap_name = base.name + str(N)
return base.from_list(cmap_name, color_list, N)
Expand Down Expand Up @@ -228,25 +231,21 @@ def CDLParse(file_name,parm):
# Allocate fortran PFT arrays
# ==============================================================================

iret=f90_pftalloc(byref(c_int(numpft)))
iret=f90_pftalloc(ci(numpft))

# ==============================================================================
# Populate the Fortran PFT structure
# ==============================================================================

# First set the arg types
f90_pftset.argtypes = \
[POINTER(c_int),POINTER(c_double),POINTER(c_int),c_char_p,c_long]

for ipft in range(numpft):
for key, parm in parms.items():
#print 'py: sending to F90: {0} = {1}'.format(parm.symbol,parm.vals[ipft])
print('{} {} '.format(parm.symbol,parm.vals[ipft]))
iret=f90_pftset(c_int(ipft+1), \
c_double(parm.vals[ipft]), \
c_int(0), \
c_char_p(parm.symbol), \
c_long(len(parm.symbol)))

c_double(parm.vals[ipft]), \
c_int(0), \
c_char_p(parm.symbol.encode('utf-8')), \
c_long(len(parm.symbol)))

# =========================================================================
# Initialize Output Arrays
Expand All @@ -259,6 +258,9 @@ def CDLParse(file_name,parm):
hd = np.zeros((numpft,ndbh))
bagwi = np.zeros((numpft,ndbh))
bagwd = np.zeros((numpft,ndbh))

bagwr = np.zeros((numpft,ndbh))

dbh = np.zeros((numpft,ndbh))
bbgw = np.zeros((numpft,ndbh))
bsapi = np.zeros((numpft,ndbh))
Expand All @@ -276,7 +278,7 @@ def CDLParse(file_name,parm):

for ipft in range(numpft):

print 'py: Solving for pft: {}'.format(ipft+1)
print('py: Solving for pft: {}'.format(ipft+1))

# Initialize Height #(d,ipft,h,dhdd)
ch_min = c_double(eparms['recruit_hgt_min'].vals[ipft])
Expand Down Expand Up @@ -316,18 +318,18 @@ def CDLParse(file_name,parm):
iret=f90_h(byref(cd),byref(cipft),byref(ch),byref(cdhdd))
hi[ipft,0] = ch.value
hd[ipft,0] = ch.value
print 'py: initialize h[{},0]={}'.format(ipft+1,ch.value)
print('py: initialize h[{},0]={}'.format(ipft+1,ch.value))

# Initialize AGB #(d,ipft,bagw,dbagwdd)
iret=f90_bagw(byref(cd),byref(cipft),byref(cbagw),byref(cdbagwdd))
bagwi[ipft,0] = cbagw.value
print 'py: initialize bagwi[{},0]={}'.format(ipft+1,cbagw.value)
print('py: initialize bagwi[{},0]={}'.format(ipft+1,cbagw.value))

# Initialize bleaf #(d,ipft,canopy_trim,bl,dbldd)
iret=f90_bleaf(byref(cd),byref(cipft),byref(ccanopy_trim),byref(cblmax),byref(cdblmaxdd))
blmaxi[ipft,0] = cblmax.value
blmaxd[ipft,0] = cblmax.value
print 'py: initialize blmaxi[{},0]={}'.format(ipft+1,cblmax.value)
print('py: initialize blmaxi[{},0]={}'.format(ipft+1,cblmax.value))

# Initialize bstore #(d,ipft,canopy_trim,bstore,dbstoredd)
iret=f90_bstore(byref(cd),byref(cipft),byref(ccanopy_trim),byref(cbstore),byref(cdbstoredd))
Expand All @@ -338,8 +340,10 @@ def CDLParse(file_name,parm):
# (dbh, nplant, site_spread, ipft, c_area,inverse)
iret= f90_carea(byref(cd),byref(cnplant),byref(csite_spread),byref(cipft),byref(ccamin),byref(cdo_reverse))
camin[ipft,0] = ccamin.value


ldense[ipft,0] = blmaxi[ipft,0]/camin[ipft,0]
print 'py: initialize careai[{},0]={}'.format(ipft+1,ccamin.value)
print('py: initialize careai[{},0]={}'.format(ipft+1,ccamin.value))

#f90_treelai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top)
cvcmax=c_double(eparms['vcmax25top'].vals[ipft])
Expand All @@ -350,30 +354,32 @@ def CDLParse(file_name,parm):
iret=f90_bfineroot(byref(cd),byref(cipft),byref(ccanopy_trim), \
byref(cbfrmax),byref(cdbfrmaxdd))
bfrmax[ipft,0] = cbfrmax.value
print 'py: initialize bfrmax[{},0]={}'.format(ipft+1,cbfrmax.value)
print('py: initialize bfrmax[{},0]={}'.format(ipft+1,cbfrmax.value))

# Initialize coarse roots #(d,ipft,bbgw,dbbgwdd)
iret=f90_bbgw(byref(cd),byref(cipft),byref(c_double(1.0)), \
byref(cbbgw),byref(cdbbgwdd))
bbgw[ipft,0] = cbbgw.value
print 'py: initialize bbgw[{},0]={}'.format(ipft+1,cbbgw.value)
print('py: initialize bbgw[{},0]={}'.format(ipft+1,cbbgw.value))


# Initialize bsap (d,ipft,canopy_trim,asapw,bsap,dbsapdd)
iret=f90_bsap(byref(cd),byref(cipft),byref(ccanopy_trim),byref(casapw),byref(cbsap),byref(cdbsapdd))
bsapi[ipft,0] = cbsap.value
bsapd[ipft,0] = cbsap.value
asapd[ipft,0] = casapw.value
print 'py: initialize bsapi[{},0]={}'.format(ipft+1,cbsap.value)
print('py: initialize bsapi[{},0]={}'.format(ipft+1,cbsap.value))

# bdead #(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeaddd)
iret=f90_bdead(byref(cbagw),byref(cbbgw),byref(cbsap),byref(cipft), \
byref(cbdead),byref(cdbagwdd),byref(cdbbgwdd), \
byref(cdbsapdd),byref(cdbdeaddd))

bdead[ipft,0] = cbdead.value
print 'py: initialize bdead[{},0]={}'.format(ipft+1,cbdead.value)
print('py: initialize bdead[{},0]={}'.format(ipft+1,cbdead.value))

bagwr[ipft,0] = (bdead[ipft,0]) * 0.6

# the metric that shan't be spoken
blmax_o_dbagwdh[ipft,0] = blmaxi[ipft,0]/(cdbagwdd.value/cdhdd.value)

Expand Down Expand Up @@ -410,7 +416,7 @@ def CDLParse(file_name,parm):
bagwi[ipft,idi] = bagwi[ipft,idi-1] + cdbagwdd.value*dd

# diagnose bleaf #(d,ipft,blmax,dblmaxdd)
iret=f90_bleaf(byref(cdc),byref(cipft),byref(c_double(1.0)),byref(cblmax),byref(cdblmaxdd))
iret=f90_bleaf(byref(cdc),byref(cipft),byref(ccanopy_trim),byref(cblmax),byref(cdblmaxdd))
blmaxd[ipft,idi] = cblmax.value

# bstore #(d,ipft,canopy_trim,bstore,dbstoredd)
Expand All @@ -426,6 +432,9 @@ def CDLParse(file_name,parm):
treelai[ipft,idi]=f90_treelai(byref(cblmax),byref(cipft),byref(ccamin), \
byref(cnplant),byref(cilayer),byref(ccanopy_lai),byref(cvcmax))




# integrate bleaf #(d,ipft,blmax,dblmaxdd)
iret=f90_bleaf(byref(cdp),byref(cipft),byref(c_double(1.0)),byref(cblmax),byref(cdblmaxdd))
blmaxi[ipft,idi] = blmaxi[ipft,idi-1] + cdblmaxdd.value*dd
Expand All @@ -450,6 +459,9 @@ def CDLParse(file_name,parm):
iret=f90_bsap(byref(cdp),byref(cipft),byref(ccanopy_trim),byref(casapw),byref(cbsap),byref(cdbsapdd))
bsapi[ipft,idi] = bsapi[ipft,idi-1] + cdbsapdd.value*dd




# the metric that shan't be spoken
# previous t-step derivatives are used for simplicity
if cdhdd.value<0.000001:
Expand All @@ -463,15 +475,17 @@ def CDLParse(file_name,parm):

# Diagnose bdead (bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeaddd)

iret=f90_bdead(byref(c_double(bagwi[ipft,idi])), \
iret=f90_bdead(byref(c_double(bagwd[ipft,idi])), \
byref(c_double(bbgw[ipft,idi])), \
byref(c_double(bsapi[ipft,idi])), \
byref(c_double(bsapd[ipft,idi])), \
byref(cipft), byref(cbdead), \
byref(cdbagwdd),byref(cdbbgwdd), \
byref(cdbsapdd),byref(cdbdeaddd))
bdead[ipft,idi] = cbdead.value


bagwr[ipft,idi] = (bdead[ipft,idi] + bsapd[ipft,idi]) * 0.6

# Create the appropriate number of line-styles, colors and widths
linestyles_base = ['-', '--', '-.', ':']
linestyles=[]
Expand Down Expand Up @@ -509,6 +523,18 @@ def CDLParse(file_name,parm):
plt.close(fig0)



if(False):
fig1_12 = plt.figure()
for ipft in range(numpft):
plt.plot(bagwd[ipft,:],bagwr[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
plt.xlabel('bagw [m]')
plt.ylabel('bagr [m]')
plt.title('')
plt.grid(True)
plt.savefig("plots/bagw_vs_bagwr.png")


if(True):
fig1 = plt.figure()
figleg = plt.figure()
Expand All @@ -520,17 +546,7 @@ def CDLParse(file_name,parm):
plt.grid(True)
plt.tight_layout()

if(False):
fig1_0 = plt.figure()
for ipft in range(numpft):
plt.plot(dbh[ipft,0:15],hi[ipft,0:15],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
plt.xlabel('diameter [cm]')
plt.ylabel('height [m]')
plt.title('Integrated Heights')
plt.grid(True)
plt.tight_layout()

if(False):
if(True):
fig1_1 = plt.figure()
for ipft in range(numpft):
plt.plot(hd[ipft,:],hi[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
Expand Down Expand Up @@ -662,7 +678,7 @@ def CDLParse(file_name,parm):
linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
ax.set_xlabel('diameter [cm]')
ax.set_ylabel('[kgC/kgC]')
ax.set_title('Sapwood (fraction of total live)')
ax.set_title('Sapwood (fraction of live)')
ax.grid(True)
# Leaf
ax = fig7_2.add_subplot(222)
Expand All @@ -671,7 +687,7 @@ def CDLParse(file_name,parm):
linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
ax.set_xlabel('diameter [cm]')
ax.set_ylabel('[kgC/kgC]')
ax.set_title('Leaf (fraction of total live)')
ax.set_title('Leaf (fraction of live)')
ax.grid(True)
# Fine Root
ax = fig7_2.add_subplot(223)
Expand All @@ -680,7 +696,7 @@ def CDLParse(file_name,parm):
linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
ax.set_xlabel('diameter [cm]')
ax.set_ylabel('[kgC/kgC]')
ax.set_title('Fine-Root (fraction of total live)')
ax.set_title('Fine-Root (fraction of live)')
ax.grid(True)
# Storage
ax = fig7_2.add_subplot(224)
Expand All @@ -689,7 +705,7 @@ def CDLParse(file_name,parm):
linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
ax.set_xlabel('diameter [cm]')
ax.set_ylabel('[kgC/kgC]')
ax.set_title('Storage (fraction of total live)')
ax.set_title('Storage (fraction of live)')
ax.grid(True)

plt.tight_layout()
Expand All @@ -698,18 +714,17 @@ def CDLParse(file_name,parm):

if(True):
fig8=plt.figure()
ax = fig8.add_subplot(111)
for ipft in range(numpft):
plt.semilogy(dbh[ipft,:],treelai[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
plt.xlabel('diameter [cm]')
plt.ylabel('[m2/m2]')
plt.title('In-Crown LAI')
plt.grid(True)
ax.plot(dbh[ipft,:],treelai[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth)
ax.ticklabel_format(style='plain')
ax.set_xlabel('diameter [cm]')
ax.set_ylabel('[m2/m2]')
ax.set_title('Untrimmed In-Crown LAI')
ax.grid(True)
plt.tight_layout()


# print(blmaxi[2,:])
# print(bfrmax[2,:])
# print(bstore[2,:])
# print(bsapd[2,:])


plt.show()
2 changes: 2 additions & 0 deletions functional_unit_testing/shared/py_src/PyF90Utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,8 @@ def ci(i8):
def cchar(fchar):
return(byref(c_char(fchar)))

def cchar3(fchar):
return(byref(c_char(fchar.encode('utf-8'))))

# We do NOT pass arrays back by reference
# This is because we will need to get their length
Expand Down