diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 3be628042b..410160ca0d 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -348,6 +348,7 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: cwd_ag_local(ncwd) ! initial value of above ground coarse woody debris. KgC/m2 real(r8) :: cwd_bg_local(ncwd) ! initial value of below ground coarse woody debris. KgC/m2 integer :: levcan ! canopy level + real(r8) :: leaf_burn_frac real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: fnrt_c ! fineroot carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] @@ -682,6 +683,33 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral nc%lmort_infra = currentCohort%lmort_infra + + ! -------------------------------------------------------------------- + ! Burn parts of trees that did *not* die in the fire. + ! currently we only remove leaves. branch and associated + ! sapwood consumption coming soon. + ! PART 4) Burn parts of grass that are consumed by the fire. + ! grasses are not killed directly by fire. They die by losing all of + ! their leaves and starving. + ! -------------------------------------------------------------------- + + leaf_c = nc%prt%GetState(leaf_organ, all_carbon_elements) + + if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then + leaf_burn_frac = currentCohort%fraction_crown_burned + else + leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) + endif + + call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) + + !KgC/gridcell/day + currentSite%flux_out = currentSite%flux_out + leaf_burn_frac* leaf_c * nc%n + + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + leaf_burn_frac * leaf_c * nc%n + + currentCohort%fraction_crown_burned = 0.0_r8 + nc%fraction_crown_burned = 0.0_r8 ! Logging is the dominant disturbance @@ -1063,9 +1091,8 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si real(r8) :: bcroot ! amount of below ground coarse root per cohort kgC. (goes into CWD_BG) real(r8) :: bstem ! amount of above ground stem biomass per cohort kgC.(goes into CWG_AG) real(r8) :: dead_tree_density ! no trees killed by fire per m2 + real(r8) :: dead_tree_num ! total number of trees killed by fire reaL(r8) :: burned_litter ! amount of each litter pool burned by fire. kgC/m2/day - real(r8) :: burned_leaves ! amount of tissue consumed by fire for leaves. KgC/individual/day - real(r8) :: leaf_burn_frac ! fraction of leaves burned real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: fnrt_c ! fineroot carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] @@ -1087,7 +1114,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si !PART 1) Burn the fractions of existing litter in the new patch that were consumed by the fire. !************************************/ do c = 1,ncwd - burned_litter = new_patch%cwd_ag(c) * patch_site_areadis/new_patch%area * & + burned_litter = currentPatch%cwd_ag(c) * patch_site_areadis/new_patch%area * & currentPatch%burnt_frac_litter(c) !kG/m2/day new_patch%cwd_ag(c) = new_patch%cwd_ag(c) - burned_litter currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day @@ -1096,7 +1123,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si enddo do p = 1,numpft - burned_litter = new_patch%leaf_litter(p) * patch_site_areadis/new_patch%area * & + burned_litter = currentPatch%leaf_litter(p) * patch_site_areadis/new_patch%area * & currentPatch%burnt_frac_litter(dl_sf) new_patch%leaf_litter(p) = new_patch%leaf_litter(p) - burned_litter currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day @@ -1129,23 +1156,28 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si bstem = (sapw_c + struct_c) * EDPftvarcon_inst%allom_agb_frac(p) ! coarse root biomass per tree bcroot = (sapw_c + struct_c) * (1.0_r8 - EDPftvarcon_inst%allom_agb_frac(p) ) - ! density of dead trees per m2. - dead_tree_density = (currentCohort%fire_mort * currentCohort%n*patch_site_areadis/currentPatch%area) / AREA + + ! Total number of dead trees + dead_tree_num = currentCohort%fire_mort * currentCohort%n*patch_site_areadis/currentPatch%area + + ! density of dead trees per m2 (spread over the new and pre-existing patch) + dead_tree_density = dead_tree_num / (new_patch%area + currentPatch%area-patch_site_areadis ) if( hlm_use_planthydro == itrue ) then - call AccumulateMortalityWaterStorage(currentSite,currentCohort,dead_tree_density*AREA) + call AccumulateMortalityWaterStorage(currentSite,currentCohort,dead_tree_num) end if ! Unburned parts of dead tree pool. ! Unburned leaves and roots - new_patch%leaf_litter(p) = new_patch%leaf_litter(p) + dead_tree_density * leaf_c * (1.0_r8-currentCohort%fraction_crown_burned) - - new_patch%root_litter(p) = new_patch%root_litter(p) + dead_tree_density * (fnrt_c+store_c) + new_patch%leaf_litter(p) = new_patch%leaf_litter(p) + & + dead_tree_density * leaf_c * (1.0_r8-currentCohort%fraction_crown_burned) currentPatch%leaf_litter(p) = currentPatch%leaf_litter(p) + dead_tree_density * & leaf_c * (1.0_r8-currentCohort%fraction_crown_burned) + new_patch%root_litter(p) = new_patch%root_litter(p) + dead_tree_density * (fnrt_c+store_c) + currentPatch%root_litter(p) = currentPatch%root_litter(p) + dead_tree_density * & (fnrt_c + store_c) @@ -1160,7 +1192,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si ! below ground coarse woody debris from burned trees do c = 1,ncwd - new_patch%cwd_bg(c) = new_patch%cwd_bg(c) + dead_tree_density * SF_val_CWD_frac(c) * bcroot + new_patch%cwd_bg(c) = new_patch%cwd_bg(c) + dead_tree_density * SF_val_CWD_frac(c) * bcroot currentPatch%cwd_bg(c) = currentPatch%cwd_bg(c) + dead_tree_density * SF_val_CWD_frac(c) * bcroot ! track as diagnostic fluxes @@ -1225,51 +1257,6 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si currentCohort => currentCohort%taller enddo ! currentCohort - - !************************************/ - ! PART 3) Burn parts of trees that did *not* die in the fire. - ! currently we only remove leaves. branch and assocaited sapwood consumption coming soon. - ! PART 4) Burn parts of grass that are consumed by the fire. - ! grasses are not killed directly by fire. They die by losing all of their leaves and starving. - !************************************/ - currentCohort => new_patch%shortest - do while(associated(currentCohort)) - - sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) - - call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area) - - if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then - burned_leaves = leaf_c * currentCohort%fraction_crown_burned - else - burned_leaves = leaf_c * currentPatch%burnt_frac_litter(lg_sf) - endif - - if (burned_leaves > 0.0_r8) then - - ! Remove burned leaves from the pool - if(leaf_c>nearzero) then - leaf_burn_frac = burned_leaves/leaf_c - else - leaf_burn_frac = 0.0_r8 - end if - call PRTBurnLosses(currentCohort%prt, leaf_organ, leaf_burn_frac) - - !KgC/gridcell/day - currentSite%flux_out = currentSite%flux_out + burned_leaves * currentCohort%n * & - patch_site_areadis/currentPatch%area * AREA - - currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm+ burned_leaves * currentCohort%n * & - patch_site_areadis/currentPatch%area * AREA - - endif - currentCohort%fraction_crown_burned = 0.0_r8 - - currentCohort => currentCohort%taller - - enddo - endif !currentPatch%fire. end subroutine fire_litter_fluxes diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 69c3c241bc..9d9de2d086 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -39,6 +39,7 @@ module FatesPlantHydraulicsMod use FatesConstantsMod, only : pi_const use FatesConstantsMod, only : cm2_per_m2 use FatesConstantsMod, only : g_per_kg + use FatesConstantsMod, only : nearzero use EDParamsMod , only : hydr_kmax_rsurf1 use EDParamsMod , only : hydr_kmax_rsurf2 @@ -2290,7 +2291,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) real(r8), parameter :: small_theta_num = 1.e-7_r8 ! avoids theta values equalling thr or ths [m3 m-3] ! hydraulics timestep adjustments for acceptable water balance error - integer :: maxiter = 1 ! maximum iterations for timestep reduction [-] + integer :: maxiter = 5 ! maximum iterations for timestep reduction [-] integer :: imult = 3 ! iteration index multiplier [-] real(r8) :: we_area_outer ! 1D plant-soil continuum water error [kgh2o m-2 individual-1] @@ -2474,13 +2475,17 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) do while(associated(ccohort)) ccohort_hydr => ccohort%co_hydr gscan_patch = gscan_patch + ccohort%g_sb_laweight - if (gscan_patch < 0._r8) then - write(fates_log(),*) 'ERROR: negative gscan_patch!' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if ccohort => ccohort%shorter enddo !cohort + ! The HLM predicted transpiration flux even though no leaves are present? + if(bc_in(s)%qflx_transp_pa(ifp) > 1.e-10_r8 .and. gscan_patchcpatch%tallest do while(associated(ccohort)) ccohort_hydr => ccohort%co_hydr @@ -2492,18 +2497,21 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ccohort_hydr%rootuptake = 0._r8 ! Relative transpiration of this cohort from the whole patch -!! qflx_rel_tran_coh = ccohort%g_sb_laweight/gscan_patch + ! [mm H2O/cohort/s] = [mm H2O / patch / s] / [cohort/patch] - qflx_tran_veg_patch_coh = bc_in(s)%qflx_transp_pa(ifp) * ccohort%g_sb_laweight/gscan_patch + if(ccohort%g_sb_laweight>nearzero) then + qflx_tran_veg_patch_coh = bc_in(s)%qflx_transp_pa(ifp) * ccohort%g_sb_laweight/gscan_patch + + qflx_tran_veg_indiv = qflx_tran_veg_patch_coh * cpatch%area* & + min(1.0_r8,cpatch%total_canopy_area/cpatch%area)/ccohort%n !AREA / ccohort%n + else + qflx_tran_veg_patch_coh = 0._r8 + qflx_tran_veg_indiv = 0._r8 + end if - qflx_tran_veg_indiv = qflx_tran_veg_patch_coh * cpatch%area* & - min(1.0_r8,cpatch%total_canopy_area/cpatch%area)/ccohort%n !AREA / ccohort%n - - ! [mm H2O/cohort/s] = [mm H2O / patch / s] / [cohort/patch] -!! qflx_tran_veg_patch_coh = qflx_trans_patch_vol * qflx_rel_tran_coh - call updateWaterDepTreeHydProps(sites(s),ccohort,bc_in(s)) - + call updateWaterDepTreeHydProps(sites(s),ccohort,bc_in(s)) + if(site_hydr%nlevsoi_hyd > 1) then ! BUCKET APPROXIMATION OF THE SOIL-ROOT HYDRAULIC GRADIENT (weighted average across layers) !call map2d_to_1d_shells(soilstate_inst, waterstate_inst, g, c, rs1(c,1), ccohort_hydr%l_aroot_layer*ccohort%n, & @@ -3471,7 +3479,7 @@ subroutine Hydraulics_1DSolve(cc_p, ft, z_node, v_node, ths_node, thr_node, kmax end if end do if(catch_nan) then - write(fates_log(),*)'EDPlantHydraulics returns nan at k = ', char(index_nan) + write(fates_log(),*)'EDPlantHydraulics returns nan at k = ', index_nan call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -3776,7 +3784,7 @@ subroutine flc_from_psi(ft, pm, psi_node, flc_node, site_hydr, bc_in ) bc_in%bsw_sisl(1), & flc_node) case default - write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = '//char(iswc) + write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = ', iswc call endrun(msg=errMsg(sourcefile, __LINE__)) end select end if @@ -3830,7 +3838,7 @@ subroutine dflcdpsi_from_psi(ft, pm, psi_node, dflcdpsi_node, site_hydr, bc_in ) bc_in%bsw_sisl(1), & dflcdpsi_node) case default - write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = '//char(iswc) + write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = ', iswc call endrun(msg=errMsg(sourcefile, __LINE__)) end select end if @@ -3878,7 +3886,7 @@ subroutine th_from_psi(ft, pm, psi_node, th_node, site_hydr, bc_in) call psi_from_th(ft, pm, th_node, psi_check ) if(psi_check > -1.e-8_r8) then - write(fates_log(),*)'bisect_pv returned positive value for water potential at pm = ', char(pm) + write(fates_log(),*)'bisect_pv returned positive value for water potential at pm = ', pm call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -3909,7 +3917,7 @@ subroutine th_from_psi(ft, pm, psi_node, th_node, site_hydr, bc_in) bc_in%watsat_sisl(1), & th_node) case default - write(fates_log(),*) 'invalid soil water characteristic function specified, iswc = '//char(iswc) + write(fates_log(),*) 'invalid soil water characteristic function specified, iswc = ', iswc call endrun(msg=errMsg(sourcefile, __LINE__)) end select end if @@ -4027,7 +4035,7 @@ subroutine psi_from_th(ft, pm, th_node, psi_node, site_hydr, bc_in) bc_in%bsw_sisl(1), & psi_node) case default - write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = '//char(iswc) + write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = ', iswc call endrun(msg=errMsg(sourcefile, __LINE__)) end select @@ -4078,7 +4086,7 @@ subroutine dpsidth_from_th(ft, pm, th_node, y, site_hydr, bc_in) bc_in%bsw_sisl(1), & y) case default - write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = '//char(iswc) + write(fates_log(),*) 'ERROR: invalid soil water characteristic function specified, iswc = ', iswc call endrun(msg=errMsg(sourcefile, __LINE__)) end select end if diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 80897f1c6c..9761f1b9e8 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -1021,8 +1021,8 @@ subroutine post_fire_mortality ( currentSite ) ! Equation 22 in Thonicke et al. 2010. currentCohort%crownfire_mort = EDPftvarcon_inst%crown_kill(currentCohort%pft)*currentCohort%fraction_crown_burned**3.0_r8 ! Equation 18 in Thonicke et al. 2010. - currentCohort%fire_mort = currentCohort%crownfire_mort+currentCohort%cambial_mort- & - (currentCohort%crownfire_mort*currentCohort%cambial_mort) !joint prob. + currentCohort%fire_mort = max(0._r8,min(1.0_r8,currentCohort%crownfire_mort+currentCohort%cambial_mort- & + (currentCohort%crownfire_mort*currentCohort%cambial_mort))) !joint prob. else currentCohort%fire_mort = 0.0_r8 !I have changed this to zero and made the mode of death removal of leaves... endif !trees diff --git a/functional_unit_testing/allometry/AutoGenVarCon.py b/functional_unit_testing/allometry/AutoGenVarCon.py new file mode 100644 index 0000000000..8879b2b7da --- /dev/null +++ b/functional_unit_testing/allometry/AutoGenVarCon.py @@ -0,0 +1,140 @@ + + +# Walk through lines of a file, if a line contains +# the string of interest (EDPftvarcon_inst), then +# parse the string to find the variable name, and save that +# to the list + + +class ParamType: + + def __init__(self,var_sym,n_dims): + + self.var_sym = var_sym + self.n_dims = n_dims + self.var_name = '' + + + + +def CheckFile(filename,check_str): + file_ptr = file(filename) + var_list = [] + found = False + for line in file_ptr: + if check_str in line: + line_split = line.split() + # substr = [i for i in line_split if check_str in i][0] + substr = line + p1 = substr.find('%')+1 + if(p1>0): + substr=substr[p1:] + p2 = substr.find('(') + p3 = substr.find(')') + # Count the number of commas between p2 and p3 + n_dims = substr[p2:p3].count(',')+1 + if(p2>0): + var_list.append(ParamType(substr[:p2],n_dims)) + + unique_list = [] + for var in var_list: + found = False + for uvar in unique_list: + if (var.var_sym == uvar.var_sym): + found = True + if(not found): + unique_list.append(var) + + return(unique_list) + + + +check_str = 'EDPftvarcon_inst%' +filename = '../../biogeochem/FatesAllometryMod.F90' + +var_list = CheckFile(filename,check_str) + + +# Add symbols here + +var_list.append(ParamType('hgt_min',1)) + + +# Now look through EDPftvarcon.F90 to determine the variable name in file +# that is associated with the variable pointer + +filename = '../../main/EDPftvarcon.F90' + +f = open(filename,"r") +contents = f.readlines() + + +var_name_list = [] +for var in var_list: + for i,line in enumerate(contents): + if (var.var_sym in line) and ('data' in line) and ('=' in line): + var.var_name = contents[i-2].split()[-1].strip('\'') + print("{} {} {}".format(var.var_sym,var.var_name,var.n_dims)) + + +f = open("f90src/AllomUnitWrap.F90_in", "r") +contents = f.readlines() +f.close() + +# Identify where we define the variables, and insert the variable definitions + +for i,str in enumerate(contents): + if 'VARIABLE-DEFINITIONS-HERE' in str: + index0=i + +index=index0+2 +for var in var_list: + if(var.n_dims==1): + contents.insert(index,' real(r8),pointer :: {}(:)\n'.format(var.var_sym)) + elif(var.n_dims==2): + contents.insert(index,' real(r8),pointer :: {}(:,:)\n'.format(var.var_sym)) + else: + print('Incorrect number of dims...') + exit(-2) + index=index+1 + +# Identify where we do the pointer assignments, and insert the pointer assignments + + +for i,str in enumerate(contents): + if 'POINTER-SPECIFICATION-HERE' in str: + index0=i + +index=index0+2 +for ivar,var in enumerate(var_list): + if(var.n_dims==1): + ins_l1='\t allocate(EDPftvarcon_inst%{}(1:numpft))\n'.format(var.var_sym) + ins_l2='\t EDPftvarcon_inst%{}(:) = nan\n'.format(var.var_sym) + ins_l3='\t iv1 = iv1 + 1\n' + ins_l4='\t EDPftvarcon_ptr%var1d(iv1)%var_name = "{}"\n'.format(var.var_name) + ins_l5='\t EDPftvarcon_ptr%var1d(iv1)%var_rp => EDPftvarcon_inst%{}\n'.format(var.var_sym) + ins_l6='\t EDPftvarcon_ptr%var1d(iv1)%vtype = 1\n' + ins_l7='\n' + if(var.n_dims==2): + ins_l1='\t allocate(EDPftvarcon_inst%{}(1:numpft,1))\n'.format(var.var_sym) + ins_l2='\t EDPftvarcon_inst%{}(:,:) = nan\n'.format(var.var_sym) + ins_l3='\t iv2 = iv2 + 1\n' + ins_l4='\t EDPftvarcon_ptr%var2d(iv2)%var_name = "{}"\n'.format(var.var_name) + ins_l5='\t EDPftvarcon_ptr%var2d(iv2)%var_rp => EDPftvarcon_inst%{}\n'.format(var.var_sym) + ins_l6='\t EDPftvarcon_ptr%var2d(iv2)%vtype = 1\n' + ins_l7='\n' + + contents.insert(index,ins_l1) + contents.insert(index+1,ins_l2) + contents.insert(index+2,ins_l3) + contents.insert(index+3,ins_l4) + contents.insert(index+4,ins_l5) + contents.insert(index+5,ins_l6) + contents.insert(index+6,ins_l7) + index=index+7 + + +f = open("f90src/AllomUnitWrap.F90", "w+") +contents = "".join(contents) +f.write(contents) +f.close() diff --git a/functional_unit_testing/allometry/drive_allomtests.py b/functional_unit_testing/allometry/drive_allomtests.py new file mode 100644 index 0000000000..9c5d5fcb38 --- /dev/null +++ b/functional_unit_testing/allometry/drive_allomtests.py @@ -0,0 +1,715 @@ +import numpy as np +import math +import matplotlib.pyplot as plt +import matplotlib as mp +import ctypes +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())) + + +# ======================================================================================= +# Set some constants. If they are used as constant arguments to the F90 routines, +# define them with their ctype identifiers +# ======================================================================================= + +ndbh = 200 +maxdbh = 50 +ccanopy_trim = c_double(1.0) # Crown Trim (0=0% of target, 1=100% of targ) +csite_spread = c_double(0.0) # Canopy spread (0=closed, 1=open) +cnplant = c_double(1.0) # Number of plants (don't change) +cilayer = c_int(1) # Index of the plant's canopy layer +ccanopy_lai = (2 * c_double)(1.0,1.0) # The LAI of the different canopy layers + # THIS VECTOR MUST MATCH ncanlayer +cdo_reverse = c_bool(0) # DO NOT GET REVERSE CROWN AREA + +# ======================================================================================= +# Setup references to fortran shared libraries +# ======================================================================================= + +allom_const_object = "./include/FatesConstantsMod.o" +allom_wrap_object = "./include/AllomUnitWrap.o" +allom_lib_object = "./include/FatesAllometryMod.o" + +# ============================================================================== +# Instantiate fortran allometry and other libraries +# ============================================================================== + +f90constlib= ctypes.CDLL(allom_const_object,mode=ctypes.RTLD_GLOBAL) +f90wraplib = ctypes.CDLL(allom_wrap_object,mode=ctypes.RTLD_GLOBAL) +f90funclib = ctypes.CDLL(allom_lib_object,mode=ctypes.RTLD_GLOBAL) + +# ======================================================================================= +# Create aliases to all of the different routines, set return types for functions +# ======================================================================================= + +f90_pftalloc = f90wraplib.__edpftvarcon_MOD_edpftvarconalloc #(numpft) +f90_pftset = f90wraplib.__edpftvarcon_MOD_edpftvarconpyset +f90_pftset.argtypes = [POINTER(c_int),POINTER(c_double),POINTER(c_int),c_char_p,c_long] +f90_h2d = f90funclib.__fatesallometrymod_MOD_h2d_allom #(h,ipft,d,dddh) +f90_h = f90funclib.__fatesallometrymod_MOD_h_allom #(d,ipft,h,dhdd) +f90_bagw = f90funclib.__fatesallometrymod_MOD_bagw_allom #(d,ipft,bagw,dbagwdd) +f90_bleaf = f90funclib.__fatesallometrymod_MOD_bleaf #(d,ipft,canopy_trim,bl,dbldd) +f90_bsap = f90funclib.__fatesallometrymod_MOD_bsap_allom #(d,ipft,canopy_trim,asapw,bsap,dbsapdd) +f90_bstore = f90funclib.__fatesallometrymod_MOD_bstore_allom #(d,ipft,canopy_trim,bstore,dbstoredd) +f90_bbgw = f90funclib.__fatesallometrymod_MOD_bbgw_allom #(d,ipft,canopy_trim,bbgw,dbbgwdd) +f90_bfineroot = f90funclib.__fatesallometrymod_MOD_bfineroot #(d,ipft,canopy_trim,bfr,dbfrdd) +f90_bdead = f90funclib.__fatesallometrymod_MOD_bdead_allom #(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeaddd) +f90_carea = f90funclib.__fatesallometrymod_MOD_carea_allom #(d,nplant,site_spread,ipft,c_area)(d,nplant,site_spread,ipft,c_area) +f90_treelai = f90funclib.__fatesallometrymod_MOD_tree_lai #(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top) +f90_treelai.restype = c_double + + +# This is the object type that holds our parameters +# ======================================================================================= +class parameter: + + def __init__(self,symbol): + + self.dtype = -9 + self.symbol = symbol + self.vals = [] + + def setval(self,val,ipft): + + self.vals[ipft] = val + +# This is just a helper script that generates random colors +# ======================================================================================= +def DiscreteCubeHelix(N): + + base = plt.cm.get_cmap('cubehelix') + np.random.seed(1) + 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) + + +# This will look through a CDL file for the provided parameter and determine +# the parameter's type, as well as fill an array with the data +# ======================================================================================= +def CDLParse(file_name,parm): + + fp = open(file_name,"r") + contents = fp.readlines() + fp.close() + + # Look in the file for the parameters + # symbol/name, record the line number + iline=-1 + isfirst = True + for i,line in enumerate(contents): + if(parm.symbol in line): + iline=i + if(isfirst): + dtype = line.split()[0] + if(dtype.strip()=="float" or (dtype.strip()=="double")): + parm.dtype = 0 + elif(dtype.strip()=="char"): + parm.dtype = 1 + isFirst=False + + if(iline==-1): + print('Could not find symbol: {} in file: {}'.format(parm.symbol,file_name)) + exit(2) + else: + search_field=True + line="" + lcount=0 + while(search_field and (lcount<100)): + line+=contents[iline] + if(line.count(';')>0): + search_field=False + else: + search_field=True + lcount=lcount+1 + iline=iline+1 + + # Parse the line + line_split = re.split(',|=',line) + # Remove the variable name entry + del line_split[0] + + # This is for read numbers + if(parm.dtype == 0): + ival=0 + for str0 in line_split: + str="" + isnum=False + for s in str0: + if (s.isdigit() or s=='.'): + str+=s + isnum=True + if(isnum): + parm.vals.append(float(str)) + + # This is a sting + elif(parm.dtype == 1): + for str0 in line_split: + # Loop several times to trim stuff off + for i in range(5): + str0=str0.strip().strip('\"').strip(';').strip() + parm.vals.append(str0) + + return(parm) + + + +# Read in the arguments +# ======================================================================================= + +parser = argparse.ArgumentParser(description='Parse command line arguments to this script.') +parser.add_argument('--fin', '--input', dest='fnamein', type=str, help="Input CDL filename. Required.", required=True) +args = parser.parse_args() + + +# Read in the parameters of interest that are used in the fortran objects. These +# parameters will be passed to the fortran allocation. +# ======================================================================================= + +parms = {} +parms['dbh_maxheight'] = CDLParse(args.fnamein,parameter('fates_allom_dbh_maxheight')) +parms['hmode'] = CDLParse(args.fnamein,parameter('fates_allom_hmode')) +parms['amode'] = CDLParse(args.fnamein,parameter('fates_allom_amode')) +parms['lmode'] = CDLParse(args.fnamein,parameter('fates_allom_lmode')) +parms['smode'] = CDLParse(args.fnamein,parameter('fates_allom_smode')) +parms['cmode'] = CDLParse(args.fnamein,parameter('fates_allom_cmode')) +parms['fmode'] = CDLParse(args.fnamein,parameter('fates_allom_fmode')) +parms['stmode'] = CDLParse(args.fnamein,parameter('fates_allom_stmode')) +parms['cushion'] = CDLParse(args.fnamein,parameter('fates_alloc_storage_cushion')) +parms['d2h1'] = CDLParse(args.fnamein,parameter('fates_allom_d2h1')) +parms['d2h2'] = CDLParse(args.fnamein,parameter('fates_allom_d2h2')) +parms['d2h3'] = CDLParse(args.fnamein,parameter('fates_allom_d2h3')) +parms['agb1'] = CDLParse(args.fnamein,parameter('fates_allom_agb1')) +parms['agb2'] = CDLParse(args.fnamein,parameter('fates_allom_agb2')) +parms['agb3'] = CDLParse(args.fnamein,parameter('fates_allom_agb3')) +parms['agb4'] = CDLParse(args.fnamein,parameter('fates_allom_agb4')) +parms['d2bl1'] = CDLParse(args.fnamein,parameter('fates_allom_d2bl1')) +parms['d2bl2'] = CDLParse(args.fnamein,parameter('fates_allom_d2bl2')) +parms['d2bl3'] = CDLParse(args.fnamein,parameter('fates_allom_d2bl3')) +parms['wood_density'] = CDLParse(args.fnamein,parameter('fates_wood_density')) +parms['c2b'] = CDLParse(args.fnamein,parameter('fates_c2b')) +parms['la_per_sa_int'] = CDLParse(args.fnamein,parameter('fates_allom_la_per_sa_int')) +parms['la_per_sa_slp'] = CDLParse(args.fnamein,parameter('fates_allom_la_per_sa_slp')) +parms['slatop'] = CDLParse(args.fnamein,parameter('fates_leaf_slatop')) +parms['slamax'] = CDLParse(args.fnamein,parameter('fates_leaf_slamax')) +parms['l2fr'] = CDLParse(args.fnamein,parameter('fates_allom_l2fr')) +parms['agb_frac'] = CDLParse(args.fnamein,parameter('fates_allom_agb_frac')) +parms['blca_expnt_diff'] = CDLParse(args.fnamein,parameter('fates_allom_blca_expnt_diff')) +parms['d2ca_coeff_min'] = CDLParse(args.fnamein,parameter('fates_allom_d2ca_coefficient_min')) +parms['d2ca_coeff_max'] = CDLParse(args.fnamein,parameter('fates_allom_d2ca_coefficient_max')) +parms['sai_scaler'] = CDLParse(args.fnamein,parameter('fates_allom_sai_scaler')) + +# Read in the parameters that are not necessary for the F90 allometry algorithms, +# but are useful for these scripts (e.g. the name of the parameter, and minimum height) +# ======================================================================================= + +eparms = {} +eparms['recruit_hgt_min'] = CDLParse(args.fnamein,parameter('fates_recruit_hgt_min')) +eparms['name'] = CDLParse(args.fnamein,parameter('fates_pftname')) +eparms['vcmax25top'] = CDLParse(args.fnamein,parameter('fates_leaf_vcmax25top')) + + +# Determine how many PFTs are here, also check to make sure that all parameters +# have the same number +# ======================================================================================= +numpft=-1 +for key, parm in parms.items(): + if( (len(parm.vals) == numpft) or (numpft==-1) ): + numpft=len(parm.vals) + else: + print('Bad length in PFT parameter') + print('parameter: {}, vals:'.format(parm.symbol),parm.vals) + + +# ============================================================================== +# Allocate fortran PFT arrays +# ============================================================================== + +iret=f90_pftalloc(byref(c_int(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]) + 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))) + + +# ========================================================================= +# Initialize Output Arrays +# ========================================================================= + +blmaxi = np.zeros((numpft,ndbh)) +blmaxd = np.zeros((numpft,ndbh)) +bfrmax = np.zeros((numpft,ndbh)) +hi = np.zeros((numpft,ndbh)) +hd = np.zeros((numpft,ndbh)) +bagwi = np.zeros((numpft,ndbh)) +bagwd = np.zeros((numpft,ndbh)) +dbh = np.zeros((numpft,ndbh)) +bbgw = np.zeros((numpft,ndbh)) +bsapi = np.zeros((numpft,ndbh)) +bsapd = np.zeros((numpft,ndbh)) +asapd = np.zeros((numpft,ndbh)) +bstore = np.zeros((numpft,ndbh)) +bdead = np.zeros((numpft,ndbh)) +dbhe = np.zeros((numpft,ndbh)) +camin = np.zeros((numpft,ndbh)) +ldense = np.zeros((numpft,ndbh)) +treelai = np.zeros((numpft,ndbh)) +blmax_o_dbagwdh = np.zeros((numpft,ndbh)) +blmax_o_dbagwdd = np.zeros((numpft,ndbh)) + + +for ipft in range(numpft): + + print 'py: Solving for pft: {}'.format(ipft+1) + + # Initialize Height #(d,ipft,h,dhdd) + ch_min = c_double(eparms['recruit_hgt_min'].vals[ipft]) + + cd = c_double(-9.0) + cdddh = c_double(-9.0) + cipft = c_int(ipft+1) + cinit = c_int(0) + + # Calculate the minimum dbh + iret=f90_h2d(byref(ch_min),byref(cipft),byref(cd),byref(cdddh)) + + # Generate a vector of diameters (use dbh) + dbh[ipft,:] = np.linspace(cd.value,maxdbh,num=ndbh) + + # Initialize various output vectors + cd = c_double(dbh[ipft,0]) + ch = c_double(-9.0) + cdhdd = c_double(-9.0) + cbagw = c_double(-9.0) + cdbagwdd = c_double(-9.0) + cblmax = c_double(-9.0) + cdblmaxdd = c_double(-9.0) + cbfrmax = c_double(-9.0) + cdbfrmaxdd = c_double(-9.0) + cbbgw = c_double(-9.0) + cdbbgwdd = c_double(-9.0) + cbsap = c_double(-9.0) + cdbsapdd = c_double(-9.0) + cbdead = c_double(-9.0) + cdbdeaddd = c_double(-9.0) + ccamin = c_double(-9.0) + casapw = c_double(-9.0) # Sapwood area + cbstore = c_double(-9.0) + cdbstoredd = c_double(-9.0) + + 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) + + # 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) + + # 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) + + # Initialize bstore #(d,ipft,canopy_trim,bstore,dbstoredd) + iret=f90_bstore(byref(cd),byref(cipft),byref(ccanopy_trim),byref(cbstore),byref(cdbstoredd)) + bstore[ipft,0] = cbstore.value + + # calculate crown area (d,nplant,site_spread,ipft,c_area) Using nplant = 1, generates units of m2 + # spread is likely 0.0, which is the value it tends towards when canopies close + # (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) + + #f90_treelai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top) + cvcmax=c_double(eparms['vcmax25top'].vals[ipft]) + treelai[ipft,0]=f90_treelai(byref(cblmax),byref(cipft),byref(ccamin), \ + byref(cnplant),byref(cilayer),byref(ccanopy_lai),byref(cvcmax)) + + # Initialize fine roots #(d,ipft,canopy_trim,bfr,dbfrdd) + 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) + + # 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) + + + # 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) + + # 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) + + # the metric that shan't be spoken + blmax_o_dbagwdh[ipft,0] = blmaxi[ipft,0]/(cdbagwdd.value/cdhdd.value) + + # the metric that shan't be spoken + blmax_o_dbagwdd[ipft,0] = blmaxi[ipft,0]/(cdbagwdd.value) + + for idi in range(1,ndbh): + + dp = dbh[ipft,idi-1] # previous position + dc = dbh[ipft,idi] # current position + dd = dc-dp + + cdp = c_double(dp) + cdc = c_double(dc) + cdbhe = c_double(-9.0) + cddedh = c_double(-9.0) + + if(ipft==2): + print("===") + + # integrate height #(d,ipft,h,dhdd) + iret=f90_h(byref(cdc),byref(cipft),byref(ch),byref(cdhdd)) + hi[ipft,idi] = hi[ipft,idi-1] + cdhdd.value*dd + + # diagnosed height + hd[ipft,idi] = ch.value + + # diagnose AGB #(d,h,ipft,bagw,dbagwdd) + iret=f90_bagw(byref(cdc),byref(cipft),byref(cbagw),byref(cdbagwdd)) + bagwd[ipft,idi] = cbagw.value + + # integrate AGB #(d,h,ipft,bagw,dbagwdd) + iret=f90_bagw(byref(cdp),byref(cipft),byref(cbagw),byref(cdbagwdd)) + 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)) + blmaxd[ipft,idi] = cblmax.value + + # bstore #(d,ipft,canopy_trim,bstore,dbstoredd) + iret=f90_bstore(byref(cdc),byref(cipft),byref(ccanopy_trim),byref(cbstore),byref(cdbstoredd)) + bstore[ipft,idi] = cbstore.value + + # calculate crown area (d,nplant,site_spread,ipft,c_area) Using nplant = 1, generates units of m2 + iret= f90_carea(byref(cdc),byref(cnplant),byref(csite_spread),byref(cipft),byref(ccamin),byref(cdo_reverse)) + camin[ipft,idi] = ccamin.value + + #f90_treelai(leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25top) + cvcmax=c_double(eparms['vcmax25top'].vals[ipft]) + 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 + + # leaf mass per square meter of crown + ldense[ipft,idi] = blmaxd[ipft,idi]/camin[ipft,idi] + + # integrate bfineroot #(d,ipft,canopy_trim,bfr,dbfrdd) + iret=f90_bfineroot(byref(cdp),byref(cipft),byref(c_double(1.0)),byref(cbfrmax),byref(cdbfrmaxdd)) + bfrmax[ipft,idi] = bfrmax[ipft,idi-1] + cdbfrmaxdd.value*dd + + # integrate bbgw #(d,h,ipft,bbgw,dbbgwdd) + iret=f90_bbgw(byref(cdp),byref(cipft),byref(cbbgw),byref(cdbbgwdd)) + bbgw[ipft,idi] = bbgw[ipft,idi-1] + cdbbgwdd.value*dd + + # diagnose bsap # (d,ipft,canopy_trim,asapw,bsap,dbsapdd) + iret=f90_bsap(byref(cdc),byref(cipft),byref(ccanopy_trim),byref(casapw),byref(cbsap),byref(cdbsapdd)) + bsapd[ipft,idi] = cbsap.value # Biomass + asapd[ipft,idi] = casapw.value # Area + + # integrate bsap + 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: + blmax_o_dbagwdh[ipft,idi] = None + else: + blmax_o_dbagwdh[ipft,idi] = blmaxi[ipft,idi-1]/(cdbagwdd.value/cdhdd.value) + + # the metric that shan't be spoken + # previous t-step derivatives are used for simplicity + blmax_o_dbagwdd[ipft,idi] = blmaxi[ipft,idi-1]/(cdbagwdd.value) + + # Diagnose bdead (bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeaddd) + + iret=f90_bdead(byref(c_double(bagwi[ipft,idi])), \ + byref(c_double(bbgw[ipft,idi])), \ + byref(c_double(bsapi[ipft,idi])), \ + byref(cipft), byref(cbdead), \ + byref(cdbagwdd),byref(cdbbgwdd), \ + byref(cdbsapdd),byref(cdbdeaddd)) + bdead[ipft,idi] = cbdead.value + + +# Create the appropriate number of line-styles, colors and widths +linestyles_base = ['-', '--', '-.', ':'] +linestyles=[] +for i in range(int(math.floor(float(numpft)/float(len(linestyles_base))))): + linestyles.extend(linestyles_base) +for i in range(numpft-len(linestyles)): + linestyles.append(linestyles_base[i]) + +my_colors = DiscreteCubeHelix(numpft) + + +mp.rcParams.update({'font.size': 14}) +mp.rcParams["savefig.directory"] = "" #os.chdir(os.path.dirname(__file__)) + +legfs = 12 +lwidth = 2.0 + +#code.interact(local=dict(globals(), **locals())) + +if(True): + fig0 = plt.figure() + figleg = plt.figure() + ax = fig0.add_subplot(111) + ax.axis("off") + ax.set_axis_off() + proxies = () + for ipft in range(numpft): + proxies = proxies + (mp.lines.Line2D([],[], \ + linestyle=linestyles[ipft], \ + color=my_colors(ipft), \ + label=eparms['name'].vals[ipft], \ + linewidth=lwidth),) + figleg.legend(handles=proxies,fontsize=12,frameon=False,labelspacing=0.25,loc='center') + plt.show(block=False) + plt.close(fig0) + + +if(True): + fig1 = plt.figure() + figleg = plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],hi[ipft,:],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): + 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): + fig1_1 = plt.figure() + for ipft in range(numpft): + plt.plot(hd[ipft,:],hi[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('height (diagnosed) [m]') + plt.ylabel('height (integrated) [m]') + plt.title('Height') + plt.grid(True) + plt.savefig("plots/hdhi.png") + +if(False): + fig2=plt.figure() + for ipft in range(numpft): + plt.plot(blmaxd[ipft,:],blmaxi[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diagnosed [kgC]') + plt.ylabel('integrated [kgC]') + plt.title('Maximum Leaf Biomass') + plt.grid(True) + plt.tight_layout() + +if(True): + fig3=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],blmaxi[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('mass [kgC]') + plt.title('Maximum Leaf Biomass') + plt.grid(True) + plt.tight_layout() + +if(True): + fig3_1=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,1:15],blmaxi[ipft,1:15],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('mass [kgC]') + plt.title('Maximum Leaf Biomass (saplings)') + plt.grid(True) + plt.tight_layout() + + +if(True): + fig4=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],camin[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('[m2] (closed canopy)') + plt.title('Crown Area') + plt.grid(True) + plt.tight_layout() + +if(True): + fig4_1=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],ldense[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('[kgC/m2] (closed canopy)') + plt.title('Leaf Mass Per Crown Area') + plt.grid(True) + plt.tight_layout() + + +if(True): + fig6=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],bagwi[ipft,:]/1000,linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('AGB [MgC]') + plt.title('Above Ground Biomass') + plt.grid(True) + plt.tight_layout() + +if(False): + fig6_1=plt.figure() + for ipft in range(numpft): + plt.plot(bagwd[ipft,:]/1000,bagwi[ipft,:]/1000,linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('AGBW deterministic [MgC]') + plt.ylabel('AGBW integrated [MgC]') + plt.title('Above Ground Biomass') + plt.grid(True) + plt.tight_layout() + +if(False): + fig5=plt.figure() + for ipft in range(numpft): + gpmask = np.isfinite(blmax_o_dbagwdh[ipft,:]) + plt.plot(dbh[ipft,gpmask],blmax_o_dbagwdh[ipft,gpmask],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('growth potential: bl/(dAGB/dh) [m]') + plt.title('Height Growth Potential') + plt.grid(True) + plt.tight_layout() + +if(False): + fig6=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],blmax_o_dbagwdd[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('diameter [cm]') + plt.ylabel('growth potential: bl/(dAGB/dd) [cm]') + plt.title('Diameter Growth Potential') + plt.grid(True) + plt.tight_layout() + +if(False): + fig7=plt.figure() + for ipft in range(numpft): + plt.plot(bsapd[ipft,:],bsapi[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('deterministic [kgC]') + plt.ylabel('integrated [kgC]') + plt.title('Sapwood Biomass') + plt.grid(True) + plt.tight_layout() + +if(False): + fig7_0=plt.figure() + for ipft in range(numpft): + plt.plot(dbh[ipft,:],bsapd[ipft,:],linestyle=linestyles[ipft],color=my_colors(ipft),linewidth=lwidth) + plt.xlabel('Diameter [cm]') + plt.ylabel('[kgC]') + plt.title('Sapwood Biomass') + plt.grid(True) + plt.tight_layout() + +if(True): + fig7_2=plt.figure(figsize=(8,6)) + # Sapwood + ax = fig7_2.add_subplot(221) + for ipft in range(numpft): + ax.plot(dbh[ipft,:],bsapd[ipft,:]/(bsapd[ipft,:]+blmaxi[ipft,:]+bfrmax[ipft,:]+bstore[ipft,:]), \ + 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.grid(True) + # Leaf + ax = fig7_2.add_subplot(222) + for ipft in range(numpft): + ax.plot(dbh[ipft,:],blmaxi[ipft,:]/(bsapd[ipft,:]+blmaxi[ipft,:]+bfrmax[ipft,:]+bstore[ipft,:]), \ + 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.grid(True) + # Fine Root + ax = fig7_2.add_subplot(223) + for ipft in range(numpft): + ax.plot(dbh[ipft,:],bfrmax[ipft,:]/(bsapd[ipft,:]+blmaxi[ipft,:]+bfrmax[ipft,:]+bstore[ipft,:]), \ + 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.grid(True) + # Storage + ax = fig7_2.add_subplot(224) + for ipft in range(numpft): + ax.plot(dbh[ipft,:],bstore[ipft,:]/(bsapd[ipft,:]+blmaxi[ipft,:]+bfrmax[ipft,:]+bstore[ipft,:]), \ + 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.grid(True) + + plt.tight_layout() + + + +if(True): + fig8=plt.figure() + 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) + plt.tight_layout() + + +# print(blmaxi[2,:]) +# print(bfrmax[2,:]) +# print(bstore[2,:]) +# print(bsapd[2,:]) + +plt.show() diff --git a/functional_unit_testing/allometry/f90src/AllomUnitWrap.F90_in b/functional_unit_testing/allometry/f90src/AllomUnitWrap.F90_in new file mode 100644 index 0000000000..471314b0bf --- /dev/null +++ b/functional_unit_testing/allometry/f90src/AllomUnitWrap.F90_in @@ -0,0 +1,218 @@ + +! ======================================================================================= +! +! This file is an alternative to key files in the fates +! filesystem. Noteably, we replace fates_r8 and fates_in +! with types that work with "ctypes". This is +! a key step in working with python +! +! We also wrap FatesGlobals to reduce the dependancy +! cascade that it pulls in from shr_log_mod. +! +! ======================================================================================= + +module shr_log_mod + + use iso_c_binding, only : c_char + use iso_c_binding, only : c_int + + contains + + function shr_log_errMsg(source, line) result(ans) + character(kind=c_char,len=*), intent(in) :: source + integer(c_int), intent(in) :: line + character(kind=c_char,len=128) :: ans + + ans = "source: " // trim(source) // " line: " + end function shr_log_errMsg + +end module shr_log_mod + + +module FatesGlobals + + contains + + integer function fates_log() + fates_log = -1 + end function fates_log + + subroutine fates_endrun(msg) + + implicit none + character(len=*), intent(in) :: msg ! string to be printed + + stop + + end subroutine fates_endrun + +end module FatesGlobals + + +module EDTypesMod + + use iso_c_binding, only : r8 => c_double + + integer, parameter :: nclmax = 2 + integer, parameter :: nlevleaf = 30 + real(r8), parameter :: dinc_ed = 1.0_r8 + +end module EDTypesMod + + +module EDPftvarcon + + use iso_c_binding, only : r8 => c_double + use iso_c_binding, only : i4 => c_int + use iso_c_binding, only : c_char + + integer,parameter :: SHR_KIND_CS = 80 ! short char + + type, public :: EDPftvarcon_inst_type + + ! VARIABLE-DEFINITIONS-HERE (DO NOT REMOVE THIS LINE, OR MOVE IT) + + end type EDPftvarcon_inst_type + + type ptr_var1 + real(r8), dimension(:), pointer :: var_rp + integer(i4), dimension(:), pointer :: var_ip + character(len=shr_kind_cs) :: var_name + integer :: vtype + end type ptr_var1 + + type ptr_var2 + real(r8), dimension(:,:), pointer :: var_rp + integer(i4), dimension(:,:), pointer :: var_ip + character(len=shr_kind_cs) :: var_name + integer :: vtype + end type ptr_var2 + + type EDPftvarcon_ptr_type + type(ptr_var1), allocatable :: var1d(:) + type(ptr_var2), allocatable :: var2d(:) + end type EDPftvarcon_ptr_type + + + type(EDPftvarcon_inst_type), public :: EDPftvarcon_inst ! ED ecophysiological constants structure + type(EDPftvarcon_ptr_type), public :: EDPftvarcon_ptr ! Pointer structure for obj-oriented id + + integer :: numparm1d ! Number of different PFT parameters + integer :: numparm2d + integer :: numpft + + logical, parameter :: debug = .true. + +contains + + + subroutine EDPftvarconPySet(ipft,rval,ival,name) + + implicit none + ! Arguments + integer(i4),intent(in) :: ipft + character(kind=c_char,len=*), intent(in) :: name + real(r8),intent(in) :: rval + integer(i4),intent(in) :: ival + ! Locals + logical :: npfound + integer :: ip + integer :: namelen + + namelen = len(trim(name)) + + if(debug) print*,"F90: ARGS: ",trim(name)," IPFT: ",ipft," RVAL: ",rval," IVAL: ",ival + + ip=0 + npfound = .true. + do ip=1,numparm1d + + if (trim(name) == trim(EDPftvarcon_ptr%var1d(ip)%var_name ) ) then + print*,"F90: Found ",trim(name)," in lookup table" + npfound = .false. + if(EDPftvarcon_ptr%var1d(ip)%vtype == 1) then ! real + EDPftvarcon_ptr%var1d(ip)%var_rp(ipft) = rval + elseif(EDPftvarcon_ptr%var1d(ip)%vtype == 2) then ! integer + EDPftvarcon_ptr%var1d(ip)%var_ip(ipft) = ival + else + print*,"F90: STRANGE TYPE" + stop + end if + end if + end do + + if(npfound)then + print*,"F90: The parameter you loaded DNE: ",name(:) + stop + end if + + do ip=1,numparm2d + if (trim(name) == trim(EDPftvarcon_ptr%var2d(ip)%var_name)) then + print*,"F90: Found ",trim(name)," in lookup table" + print*,"BUT... WE AVOID USING 2D VARIABLES FOR NOW..." + print*,"REMOVE THIS TEST" + stop + end if + end do + + + ! Perform a check to see if the target array is being filled + if (trim(name) == 'fates_allom_d2h1') then + if (EDPftvarcon_inst%allom_d2h1(ipft) == rval) then + print*,"F90: POINTER CHECK PASSES:",rval," = ",EDPftvarcon_inst%allom_d2h1(ipft) + else + print*,"F90: POINTER CHECK FAILS:",rval," != ",EDPftvarcon_inst%allom_d2h1(ipft) + stop + end if + end if + + if (trim(name) == 'fates_wood_density' ) then + if (EDPftvarcon_inst%wood_density(ipft) == rval) then + print*,"F90: POINTER CHECK PASSES:",rval," = ",EDPftvarcon_inst%wood_density(ipft) + else + print*,"F90: POINTER CHECK FAILS:",rval," != ",EDPftvarcon_inst%wood_density(ipft) + stop + end if + end if + + return + end subroutine EDPftvarconPySet + + + subroutine EDPftvarconAlloc(numpft_in) + ! + + ! !ARGUMENTS: + integer(i4), intent(in) :: numpft_in + ! LOCALS: + integer :: iv1 ! The parameter incrementer + integer :: iv2 + !------------------------------------------------------------------------ + + numpft = numpft_in + + allocate( EDPftvarcon_ptr%var1d(100)) ! Make this plenty large + allocate( EDPftvarcon_ptr%var2d(100)) + iv1=0 + iv2=0 + + ! POINTER-SPECIFICATION-HERE (DO NOT REMOVE THIS LINE, OR MOVE IT) + +! allocate( EDPftvarcon_inst%allom_dbh_maxheight (1:numpft)); EDPftvarcon_inst%allom_dbh_maxheight (:) = nan +! iv = iv + 1 +! EDPftvarcon_ptr%var1d(iv)%var_name = "fates_allom_dbh_maxheight" +! EDPftvarcon_ptr%var1d(iv)%var_rp => EDPftvarcon_inst%allom_dbh_maxheight +! EDPftvarcon_ptr%var1d(iv)%vtype = 1 + + + numparm1d = iv1 + numparm2d = iv2 + + + print*,"F90: ALLOCATED ",numparm1d," PARAMETERS, FOR ",numpft," PFTs" + + + return + end subroutine EDPftvarconAlloc + +end module EDPftvarcon diff --git a/functional_unit_testing/allometry/include/README b/functional_unit_testing/allometry/include/README new file mode 100644 index 0000000000..bfa612f78d --- /dev/null +++ b/functional_unit_testing/allometry/include/README @@ -0,0 +1 @@ +This holds the place of the include folder \ No newline at end of file diff --git a/functional_unit_testing/allometry/plots/README b/functional_unit_testing/allometry/plots/README new file mode 100644 index 0000000000..c32df9df9a --- /dev/null +++ b/functional_unit_testing/allometry/plots/README @@ -0,0 +1 @@ +Placeholder for the folder \ No newline at end of file diff --git a/functional_unit_testing/allometry/simple_build.sh b/functional_unit_testing/allometry/simple_build.sh new file mode 100755 index 0000000000..a06b9724b0 --- /dev/null +++ b/functional_unit_testing/allometry/simple_build.sh @@ -0,0 +1,57 @@ +#!/bin/bash + +FC='gfortran -g -shared -fPIC' + +# First copy over the FatesConstants file, but change the types of the fates_r8 and fates_int + +old_fates_r8_str=`grep -e integer ../../main/FatesConstantsMod.F90 | grep fates_r8 | sed 's/^[ \t]*//;s/[ \t]*$//'` +new_fates_r8_str='use iso_c_binding, only: fates_r8 => c_double' + +old_fates_int_str=`grep -e integer ../../main/FatesConstantsMod.F90 | grep fates_int | sed 's/^[ \t]*//;s/[ \t]*$//'` +new_fates_int_str='use iso_c_binding, only: fates_int => c_int' + +# Add the new lines (need position change, don't swap) + +sed "/implicit none/i $new_fates_r8_str" ../../main/FatesConstantsMod.F90 > f90src/FatesConstantsMod.F90 +sed -i "/implicit none/i $new_fates_int_str" f90src/FatesConstantsMod.F90 + +# Delete the old lines + +sed -i "/$old_fates_r8_str/d" f90src/FatesConstantsMod.F90 +sed -i "/$old_fates_int_str/d" f90src/FatesConstantsMod.F90 + + +# This re-writes the wrapper so that it uses all the correct parameters +# in FatesAllometryMod.F90 +python AutoGenVarCon.py + + +# Procedure for auto-generating AllomUnitWrap +# 1) scan FatesAllometry and create list of EDPftVarcon_inst variables +# 2) scan EDpftVarcon and get the name of the in-file parameter names associated +# with these variables + + + + +rm -f include/*.o +rm -f include/*.mod + + +# Build the new file with constants + +${FC} -I include/ -J include/ -o include/FatesConstantsMod.o f90src/FatesConstantsMod.F90 + +${FC} -I include/ -J include/ -o include/AllomUnitWrap.o f90src/AllomUnitWrap.F90 + +${FC} -I include/ -J include/ -o include/FatesAllometryMod.o ../../biogeochem/FatesAllometryMod.F90 + + +#${FC} -g -o include/FatesConstantsMod.o ../main/FatesConstantsMod.F90 + +#gfortran -shared -fPIC -g -o include/EDTypesMod.o ../main/EDTypesMod.F90 + + + + +#gfortran diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index e25987e576..e10944e898 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -2,14 +2,14 @@ module EDParamsMod ! ! module that deals with reading the ED parameter file ! - + + use FatesConstantsMod, only : r8 => fates_r8 use FatesParametersInterface, only : param_string_length use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg - use shr_kind_mod , only: r8 => shr_kind_r8 implicit none private diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 432cd85489..dd77c2cb20 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -4,7 +4,6 @@ module EDTypesMod use FatesConstantsMod, only : ifalse use FatesConstantsMod, only : itrue use FatesGlobals, only : fates_log - use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use FatesHydraulicsMemMod, only : ed_cohort_hydr_type use FatesHydraulicsMemMod, only : ed_site_hydr_type use PRTGenericMod, only : prt_vartypes diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 6c790a22b4..fa714aedcb 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -73,6 +73,9 @@ module FatesInventoryInitMod ! defined in model memory and a physical ! site listed in the file + logical, parameter :: do_inventory_out = .false. + + public :: initialize_sites_by_inventory contains @@ -127,6 +130,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) integer, allocatable :: inv_format_list(:) ! list of format specs character(len=path_strlen), allocatable :: inv_css_list(:) ! list of css file names character(len=path_strlen), allocatable :: inv_pss_list(:) ! list of pss file names + real(r8), allocatable :: inv_lat_list(:) ! list of lat coords real(r8), allocatable :: inv_lon_list(:) ! list of lon coords integer :: invsite ! index of inventory site @@ -138,13 +142,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) real(r8) :: basal_area_postf ! basal area before fusion (m2/ha) real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) - real(r8), parameter :: max_ba_diff = 1.0e-2 ! 1% is the maximum allowable - ! change in BA due to fusion - ! I. Load the inventory list file, do some file handle checks ! ------------------------------------------------------------------------------------------ sitelist_file_unit = shr_file_getUnit() + + inquire(file=trim(hlm_inventory_ctrl_file),exist=lexist,opened=lopen) if( .not.lexist ) then ! The inventory file list DNE write(fates_log(), *) 'An inventory Initialization was requested.' @@ -517,18 +520,14 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) write(fates_log(),*) basal_area_postf,' [m2/ha]' write(fates_log(),*) '-------------------------------------------------------' - ! Check to see if the fusion process has changed too much - ! We are sensitive to fusion in inventories because we may be asking for a massive amount - ! of fusion. For instance some init files are directly from inventory, where a cohort - ! is synomomous with a single plant. - - if( abs(basal_area_postf-basal_area_pref)/basal_area_pref > max_ba_diff ) then - write(fates_log(),*) 'Inventory Fusion Changed total biomass beyond reasonable limit' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end do + ! If this is flagged as true, the post-fusion inventory will be written to file + ! in the run directory. + if(do_inventory_out)then + call write_inventory_type1(sites(s)) + end if + end do deallocate(inv_format_list, inv_pss_list, inv_css_list, inv_lat_list, inv_lon_list) return @@ -1053,4 +1052,101 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & return end subroutine set_inventory_edcohort_type1 + ! ==================================================================================== + + subroutine write_inventory_type1(currentSite) + + ! -------------------------------------------------------------------------------- + ! This subroutine writes the cohort/patch inventory type files in the "type 1" + ! format. Note that for compatibility with ED2, we chose an old type that has + ! both extra unused fields and is missing fields from FATES. THis is not + ! a recommended file type for restarting a run. + ! The files will have a lat/long tag added to their name, and will be + ! generated in the run folder. + ! -------------------------------------------------------------------------------- + + use shr_file_mod, only : shr_file_getUnit + use shr_file_mod, only : shr_file_freeUnit + + ! Arguments + type(ed_site_type), target :: currentSite + + ! Locals + type(ed_patch_type), pointer :: currentpatch + type(ed_cohort_type), pointer :: currentcohort + + character(len=128) :: pss_name_out ! output file string + character(len=128) :: css_name_out ! output file string + integer :: pss_file_out + integer :: css_file_out + integer :: ilat_int,ilat_dec ! for output string parsing + integer :: ilon_int,ilon_dec ! for output string parsing + character(len=32) :: patch_str + character(len=32) :: cohort_str + integer :: ipatch + integer :: icohort + character(len=1) :: ilat_sign,ilon_sign + + ! Generate pss/css file name based on the location of the site + ilat_int = abs(int(currentSite%lat)) + ilat_dec = int(100000*(abs(currentSite%lat) - real(ilat_int,r8))) + ilon_int = abs(int(currentSite%lon)) + ilon_dec = int(100000*(abs(currentSite%lon) - real(ilon_int,r8))) + + if(currentSite%lat>=0._r8)then + ilat_sign = 'N' + else + ilat_sign = 'S' + end if + if(currentSite%lon>=0._r8)then + ilon_sign = 'E' + else + ilon_sign = 'W' + end if + + write(pss_name_out,'(A8,I2.2,A1,I5.5,A1,A1,I3.3,A1,I5.5,A1,A4)') & + 'pss_out_',ilat_int,'.',ilat_dec,ilat_sign,'_',ilon_int,'.',ilon_dec,ilon_sign,'.txt' + write(css_name_out,'(A8,I2.2,A1,I5.5,A1,A1,I3.3,A1,I5.5,A1,A4)') & + 'css_out_',ilat_int,'.',ilat_dec,ilat_sign,'_',ilon_int,'.',ilon_dec,ilon_sign,'.txt' + + pss_file_out = shr_file_getUnit() + css_file_out = shr_file_getUnit() + + open(unit=pss_file_out,file=trim(pss_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') + open(unit=css_file_out,file=trim(css_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') + + write(pss_file_out,*) 'time patch trk age area water fsc stsc stsl ssc psc msn fsn' + write(css_file_out,*) 'time patch cohort dbh hite pft nplant bdead alive Avgrg' + + ipatch=0 + currentpatch => currentSite%youngest_patch + do while(associated(currentpatch)) + ipatch=ipatch+1 + + write(patch_str,'(A7,i4.4,A)') '' + + write(pss_file_out,*) '0000 ',trim(patch_str),' 2 ',currentPatch%age,currentPatch%area/AREA, & + '0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000' + + icohort=0 + currentcohort => currentpatch%tallest + do while(associated(currentcohort)) + icohort=icohort+1 + write(cohort_str,'(A7,i4.4,A)') '' + write(css_file_out,*) '0000 ',trim(patch_str),' ',trim(cohort_str), & + currentCohort%dbh,0.0,currentCohort%pft,currentCohort%n/currentPatch%area,0.0,0.0,0.0 + + currentcohort => currentcohort%shorter + end do + currentPatch => currentpatch%older + enddo + + close(css_file_out) + close(pss_file_out) + + call shr_file_freeUnit(css_file_out) + call shr_file_freeUnit(pss_file_out) + + end subroutine write_inventory_type1 + end module FatesInventoryInitMod