Skip to content

Commit

Permalink
general improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
@rgknox
rgknox committed Oct 5, 2016
1 parent 78c605c commit 17a2064
Show file tree
Hide file tree
Showing 3 changed files with 93 additions and 41 deletions.
60 changes: 35 additions & 25 deletions EDAllomMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -806,9 +806,9 @@ subroutine d2blmax_2pwr(d,p1,p2,c2b,blmax,dblmaxdd)
!
! bl = exp(a2) * ((1/exp(a1)) * d**(1/b1))**b2
! bl = exp(a2) * (1/exp(a1))**b2 * d**(b2/b1)
! bl = p1 * d ** d2bl2_ad
! where: p1 = exp(a2) * (1/exp(a1))**b2
! d2bl2_ad = (b2/b1)
! bl = p1 * d ** p2
! where: p1 = exp(a2) * (1/exp(a1))**b2 = 0.0201198
! p2 = (b2/b1) = 3.1791044
! For T. tuberculata (canopy tree):
! a1 = -0.0704, b1 = 0.67
! a2 = -4.056, b2 = 2.13
Expand Down Expand Up @@ -862,21 +862,34 @@ subroutine dh2blmax_2pwr(d,ipft,p1,p2,c2b,blmax,dblmaxdd)
real(r8) :: dddh_eff ! effective diameter to height differential
real(r8) :: ddeffdd ! effective diameter to diameter differential

! This call is needed to calculate the rate of change of
! the actual h with d
call h_allom(d,ipft,h,dhdd)

! This call is needed to calculate the effective dbh
call h2d_allom(h,ipft,dbh_eff,dddh_eff)
associate( &
d_adult => EDecophyscon%d_adult(ipft))

blmax = p1*dbh_eff**p2/c2b

! This is the rate of change of the effective diameter
! with respect to the actual diameter (1.0 in non-height capped)
ddeffdd = dddh_eff * dhdd

dblmaxdd = p1*p2*dbh_eff**(p2-1.0_r8) / c2b * ddeffdd



! This call is needed to calculate the effective dbh
! Note that this is only necessary for adult plants
! Because only adult plants are nearing their asymptotic heights
if (d>=d_adult) then

! This call is needed to calculate the rate of change of
! the actual h with d
call h_allom(d,ipft,h,dhdd)

call h2d_allom(h,ipft,dbh_eff,dddh_eff)

! This is the rate of change of the effective diameter
! with respect to the actual diameter (1.0 in non-height capped)
ddeffdd = dddh_eff * dhdd
else
! In this case we assume that dbh_eff == d
! because we are not in some forced asymptotic portion of the curve
dbh_eff = d
ddeffdd = 1.0_r8
end if
blmax = p1*dbh_eff**p2/c2b
dblmaxdd = p1*p2*dbh_eff**(p2-1.0_r8) / c2b * ddeffdd
end associate
return
end subroutine dh2blmax_2pwr

Expand Down Expand Up @@ -940,8 +953,6 @@ subroutine d2h_chave2014(d,p1,p2,p3,eclim,dbh_maxh,h,dhdd)
! integral incredibly messy. Thus we use an Euler step, yes ugly,
! but it is a simple function so don't over-think it

print*,'STARTING CHAVE2014',d,dbh_maxh

if (d>0.5_r8*dbh_maxh) then
dbh0=0.5_r8*dbh_maxh
h = exp( p1 - eclim + p2*log(dbh0) + p3*log(dbh0)**2.0_r8 )
Expand All @@ -962,7 +973,6 @@ subroutine d2h_chave2014(d,p1,p2,p3,eclim,dbh_maxh,h,dhdd)
! display("this is innefficient, and was not thought to had been")
! display("necessary for production runs")
else
print*,"logd:",d,log(d)
h = exp( p1 - eclim + p2*log(d) + p3*log(d)**2.0 )
end if

Expand Down Expand Up @@ -990,7 +1000,6 @@ subroutine d2h_chave2014(d,p1,p2,p3,eclim,dbh_maxh,h,dhdd)
exp(p3*log(d)**2.0_r8) )
dhdd = dhpdd*(1.0_r8-fl)

print*,'END CHAVE2014'
return

end subroutine d2h_chave2014
Expand Down Expand Up @@ -1060,7 +1069,7 @@ subroutine d2h_2pwr(d,p1,p2,dbh_maxh,h,dhdd)
! log(d) = a + b*log(h)
! d = exp(a) * h**b
! h = (1/exp(a)) * d**(1/b)
! h = p1*d**p2 where p1 = 1/exp(a) p2 = 1/b
! h = p1*d**p2 where p1 = 1/exp(a) = 1.07293 p2 = 1/b = 1.4925
! d = (h/p1)**(1/p2)
! For T. tuberculata (canopy tree) a = -0.0704, b = 0.67
! =========================================================================
Expand Down Expand Up @@ -1326,11 +1335,12 @@ subroutine dh2bag_salda(d,h,ipft,bag,dbagdd)
real(r8),intent(in) :: d ! plant diameter [cm]
real(r8),intent(in) :: h ! plant height [m]
integer(li),intent(in) :: ipft ! PFT index
real(r8),intent(out) :: bag ! plant height [m]
real(r8),intent(out) :: bag ! plant biomass [kgC/indiv]
real(r8),intent(out) :: dbagdd ! change in agb per diameter [kgC/cm]

real(r8) :: term1,term2,term3,hj,dhdd

real(r8),parameter :: tot2above = 0.7 ! convert total to above fraction
real(r8),parameter :: d2bag1 = 0.06896_r8
real(r8),parameter :: d2bag2 = 0.572_r8
real(r8),parameter :: d2bag3 = 1.94_r8
Expand All @@ -1340,14 +1350,14 @@ subroutine dh2bag_salda(d,h,ipft,bag,dbagdd)
c2b => EDecophyscon%c2b(ipft), &
wood_density => EDecophyscon%wood_density(ipft))

bag = d2bag1*(h**d2bag2)*(d**d2bag3)*(wood_density**d2bag4)/c2b
bag = tot2above*d2bag1*(h**d2bag2)*(d**d2bag3)*(wood_density**d2bag4)

! bag = a1 * h**a2 * d**a3 * r**a4
! dbag/dd = a1*r**a4 * d/dd (h**a2*d**a3)
! dbag/dd = a1*r**a4 * [ d**a3 *d/dd(h**a2) + h**a2*d/dd(d**a3) ]
! dbag/dd = a1*r**a4 * [ d**a3 * a2*h**(a2-1)dh/dd + h**a2*a3*d**(a3-1)]

term1 = d2bag1*(wood_density**d2bag4)/c2b
term1 = d2bag1*(wood_density**d2bag4)
term2 = (h**d2bag2)*d2bag3*d**(d2bag3-1.0_r8)

call h_allom(d,ipft,hj,dhdd)
Expand Down
47 changes: 40 additions & 7 deletions allom_params.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,16 @@
<eclim> 0.0 </eclim>
<llspan> 1.0 </llspan>
<bl_min> 0.4 </bl_min>
<h_max> 35.0 </h_max>
<h_min> 1.25 </h_min>
<h_min> 1.00 </h_min>
<slatop> 0.024 </slatop>
<d_adult> 10.0 </d_adult>
<d_sap> 5.0 </d_sap>
<d_sap> 2.0 </d_sap>
<f2l_ratio> 1.0 </f2l_ratio>
<agb_fraction> 0.65 </agb_fraction>

<pfts>

<pft tag="martinezcano-lescure">
<hallom_mode> 5 </hallom_mode>
<hallom_ad_mode> 5 </hallom_ad_mode>
<hallom_sap_mode> 5 </hallom_sap_mode>
<lallom_ad_mode> 2 </lallom_ad_mode>
Expand All @@ -25,14 +23,15 @@
<aallom_mode> 1 </aallom_mode>
<callom_mode> 1 </callom_mode>
<sallom_mode> 1 </sallom_mode>
<h_max> 35.0 </h_max>
<latosa_int> 0.75 </latosa_int>
<latosa_slp> 0.0 </latosa_slp>
<d2h1_ad> 56.8 </d2h1_ad>
<d2h2_ad> 0.74 </d2h2_ad>
<d2h3_ad> 21.7 </d2h3_ad>
<d2h1_sap> 55.0 </d2h1_sap>
<d2h2_sap> 0.7 </d2h2_sap>
<d2h3_sap> 21.0 </d2h3_sap>
<d2h1_sap> 56.8 </d2h1_sap>
<d2h2_sap> 0.74 </d2h2_sap>
<d2h3_sap> 21.7 </d2h3_sap>
<d2bl1_ad> 0.00873</d2bl1_ad>
<d2bl2_ad> 2.1360 </d2bl2_ad>
<d2bl3_ad> -999 </d2bl3_ad>
Expand All @@ -44,6 +43,36 @@
<wood_density> 0.7 </wood_density>
</pft>

<pft tag="martinezcano-king-lescure">
<hallom_ad_mode> 5 </hallom_ad_mode>
<hallom_sap_mode> 3 </hallom_sap_mode>
<lallom_ad_mode> 2 </lallom_ad_mode>
<lallom_sap_mode> 2 </lallom_sap_mode>
<fallom_mode> 1 </fallom_mode>
<aallom_mode> 1 </aallom_mode>
<callom_mode> 1 </callom_mode>
<sallom_mode> 1 </sallom_mode>
<h_max> 35.0 </h_max>
<latosa_int> 0.75 </latosa_int>
<latosa_slp> 0.0 </latosa_slp>
<d2h1_ad> 56.8 </d2h1_ad>
<d2h2_ad> 0.74 </d2h2_ad>
<d2h3_ad> 21.7 </d2h3_ad>
<d2h1_sap> 1.07293 </d2h1_sap>
<d2h2_sap> 1.4925 </d2h2_sap>
<d2h3_sap> -999 </d2h3_sap>
<d2bl1_ad> 0.00873</d2bl1_ad>
<d2bl2_ad> 2.1360 </d2bl2_ad>
<d2bl3_ad> -999 </d2bl3_ad>
<d2bl1_sap> 0.0201 </d2bl1_sap>
<d2bl2_sap> 3.179 </d2bl2_sap>
<d2bl3_sap> -999 </d2bl3_sap>
<d2bag1> 0.0673 </d2bag1>
<d2bag2> 0.976 </d2bag2>
<wood_density> 0.7 </wood_density>
</pft>


<pft tag="chave-lescure">
<hallom_ad_mode> 1 </hallom_ad_mode>
<hallom_sap_mode> 1 </hallom_sap_mode>
Expand All @@ -53,6 +82,7 @@
<aallom_mode> 1 </aallom_mode>
<callom_mode> 1 </callom_mode>
<sallom_mode> 1 </sallom_mode>
<h_max> 35.0 </h_max>
<latosa_int> 0.75 </latosa_int>
<latosa_slp> 0.0 </latosa_slp>
<d2h1_ad> 0.893 </d2h1_ad>
Expand Down Expand Up @@ -81,6 +111,7 @@
<aallom_mode> 1 </aallom_mode>
<callom_mode> 1 </callom_mode>
<sallom_mode> 1 </sallom_mode>
<h_max> 35.0 </h_max>
<latosa_int> 0.75 </latosa_int>
<latosa_slp> 0.0 </latosa_slp>
<d2h1_ad> 0.893 </d2h1_ad>
Expand Down Expand Up @@ -111,6 +142,7 @@
<aallom_mode> 2 </aallom_mode>
<callom_mode> 1 </callom_mode>
<sallom_mode> 1 </sallom_mode>
<h_max> 35.0 </h_max>
<latosa_int> 0.75 </latosa_int>
<latosa_slp> 0.0 </latosa_slp>
<d2h1_ad> 1.7882 </d2h1_ad>
Expand Down Expand Up @@ -139,6 +171,7 @@
<aallom_mode> 3 </aallom_mode>
<callom_mode> 1 </callom_mode>
<sallom_mode> 1 </sallom_mode>
<h_max> 35.0 </h_max>
<latosa_int> 0.75 </latosa_int>
<latosa_slp> 0.0 </latosa_slp>
<d2h1_ad> 0.64 </d2h1_ad>
Expand Down
27 changes: 18 additions & 9 deletions drive_allomtests.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
from ctypes import * #byref, cdll, c_int, c_double, c_char_p, c_long
import xml.etree.ElementTree as ET

pft_xml_file = "allom_params.xml"
pft_xml_file = "params1.xml"
allom_wrap_object = "./EDAllomUnitWrap.o"
allom_lib_object = "./EDAllomMod.o"

Expand Down Expand Up @@ -409,21 +409,30 @@ def wait():
plt.savefig("plots/blmaxi.png")


fig4=plt.figure()
ax1 = fig4.add_subplot(1,2,1)
ax2 = fig4.add_subplot(1,2,2)
ax2.set_yscale('log')
#fig4=plt.figure()
#ax1 = fig4.add_subplot(1,2,1)
#ax2 = fig4.add_subplot(1,2,2)
#ax2.set_yscale('log')
#for ipft in range(numpft):
# ax1.plot(dbh[ipft,:],bagi[ipft,:],label="{}".format(pftparms[ipft]['name']))
# ax2.plot(dbh[ipft,:],bagi[ipft,:],label="{}".format(pftparms[ipft]['name']))
#plt.legend(loc='upper left')
#plt.xlabel('diameter [cm]')
#plt.ylabel('mass [kgC]')
#plt.title('Above Ground Biomass')
#plt.grid(True)
#plt.savefig("plots/agbi.png")

fig6=plt.figure()
for ipft in range(numpft):
ax1.plot(dbh[ipft,:],bagi[ipft,:],label="{}".format(pftparms[ipft]['name']))
ax2.plot(dbh[ipft,:],bagi[ipft,:],label="{}".format(pftparms[ipft]['name']))
plt.plot(dbh[ipft,:],bagi[ipft,:]/1000,label="{}".format(pftparms[ipft]['name']))
plt.legend(loc='upper left')
plt.xlabel('diameter [cm]')
plt.ylabel('mass [kgC]')
plt.ylabel('AGB [MgC]')
plt.title('Above Ground Biomass')
plt.grid(True)
plt.savefig("plots/agbi.png")


fig5=plt.figure()
for ipft in range(numpft):
gpmask = np.isfinite(blmax_o_dbagdh[ipft,:])
Expand Down

0 comments on commit 17a2064

Please sign in to comment.