Skip to content

Commit

Permalink
Merge pull request #73 from noaa-oar-arl/bio_opts
Browse files Browse the repository at this point in the history 
Bio opts
  • Loading branch information
@drnimbusrain
drnimbusrain authored Jun 12, 2023
2 parents 80c84d7 + f851b4e commit 1a5dabc
Show file tree
Hide file tree
Showing 8 changed files with 502 additions and 78 deletions.
7 changes: 7 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,6 +197,13 @@ https://nacc-in-the-cloud.s3.amazonaws.com/inputs/geo-files/gfs.canopy.t12z.2022
| `lu_opt` | integer for input model land use type (`0`: VIIRS 17 Cat (default) or `1`: MODIS-IGBP 20 Cat (valid LU types 1-10 and 12); input mapped to Massman et al.) |
| `z0_opt` | integer (`0`: use model input or `1`: vegtype dependent z0 for first estimate) |
| `bio_cce` | user-set real value of MEGAN biogenic emissions "canopy environment coefficient" used to tune emissions to model inputs/calculations (default: `0.21`, based on Silva et al. 2020) |
| `biovert_opt` | user set biogenic vertical summing option (`0`: no sum, full leaf-level biogenic emissions, units=kg/m3/s; `1`: MEGANv3-like summing of LAD weighted activity coefficients using the canopy-app plant distributions, caution-- units=kg m-2 s-1 and puts the total emissions in the topmost canopy-app model layer only; `2`: Same as in option 1, but instead uses Gaussian/normally weighted activity coefficients acoss all sub-canopy layers -- also units of kg m-2 s-1 in topmost model layer; `3`: Same as in option 1, but instead uses evenly weighted activity coefficients acoss all sub-canopy layers -- also units of kg m-2 s-1 in topmost model layer |
| `ssg_opt` | integer for using either input data (= `0`, default) or user set shrub/savanna/grass (SSG) vegetation type heights from namelist (= `1`). Currently, GEDI FCH input data only provides canopy heights for forests and not SSG. Warning: use of ssg_opt=1 will overide typically higher resolution input data (e.g., GEDI) forest canopy heights where the lower resolution vegtype data indicates SSG |
| `ssg_set` | user-set real value of constant SSG vegetation type heights (m) (only used if `ssg_opt=1`) |
| `crop_opt` | integer for using either input data (= `0`, default) or user set crop vegetation type heights from namelist (= `1`). Currently, GEDI FCH input data only provides canopy heights for forests and not crops. Warning: use of crop_opt=1 will overide typically higher resolution input data (e.g., GEDI) forest canopy heights where the lower resolution vegtype data indicates crops |
| `crop_set` | user-set real value of constant crop vegetation type heights (m) (only used if `crop_opt=1`) |
| `co2_opt` | user-set options for applying a CO2 inhibition factor for biogenic isoprene-only emissions using either the [Possell & Hewitt (2011)](https://doi.org/10.1111/j.1365-2486.2010.02306.x) (= `0`, default) or [Wilkinson et al. (2009)](https://doi.org/10.1111/j.1365-2486.2008.01803.x) method (= `1`). Use of option = `1` (Possell & Hewitt 2011) is especially recommended for sub-ambient CO2 concentrations. To turn off co2 inhibition factor set `co2_opt=2` |
| `co2_set` | user-set real value of atmospheric co2 concentration (ppmv) (only used if `co2_opt=0` or `co2_opt=1`) |
| `lai_thresh` | user-set real value of LAI threshold for contiguous canopy (m2/m2) |
| `frt_thresh` | user-set real value of forest fraction threshold for contiguous canopy |
| `fch_thresh` | user-set real value of canopy height threshold for contiguous canopy (m) |
Expand Down
63 changes: 35 additions & 28 deletions input/namelist.canopy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,32 +5,39 @@
/

&USERDEFS
infmt_opt = 0
nlat = 43
nlon = 86
modlays = 100
modres = 0.5
ifcanwind = .TRUE.
ifcanwaf = .TRUE.
ifcaneddy = .TRUE.
ifcanphot = .TRUE.
ifcanbio = .TRUE.
href_opt = 0
href_set = 10.0
z0ghc = 0.0025
rsl_opt = 0
lambdars = 1.25
dx_opt = 0
dx_set = 12000.0
flameh_opt = 2
flameh_set = 1.0
frp_fac = 1.0
pai_opt = 0
pai_set = 4.0
lu_opt = 1
z0_opt = 0
bio_cce = 0.21
lai_thresh = 0.1
frt_thresh = 0.1
fch_thresh = 0.5
infmt_opt = 0
nlat = 43
nlon = 86
modlays = 100
modres = 0.5
ifcanwind = .TRUE.
ifcanwaf = .TRUE.
ifcaneddy = .TRUE.
ifcanphot = .TRUE.
ifcanbio = .TRUE.
href_opt = 0
href_set = 10.0
z0ghc = 0.0025
rsl_opt = 0
lambdars = 1.25
dx_opt = 0
dx_set = 12000.0
flameh_opt = 2
flameh_set = 1.0
frp_fac = 1.0
pai_opt = 0
pai_set = 4.0
lu_opt = 1
z0_opt = 0
bio_cce = 0.21
biovert_opt = 0
ssg_opt = 0
ssg_set = 1.0
crop_opt = 0
crop_set = 3.0
co2_opt = 0
co2_set = 400.0
lai_thresh = 0.1
frt_thresh = 0.1
fch_thresh = 0.5
/
58 changes: 58 additions & 0 deletions python/examples.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -760,6 +760,64 @@
" ax.plot(clu, ds_c[vn])\n",
" ax.set(ylabel=vn, xlabel=\"CLU\")"
]
},
{
"cell_type": "markdown",
"id": "13cc57aa-666c-4f27-ac01-400838ba94a0",
"metadata": {},
"source": [
"### Emissions calculation options"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "578f31ab-832e-4303-88c4-9cb444fe99f9",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"cases = config_cases(\n",
" file_vars=\"../input/input_variables_point.txt\",\n",
" infmt_opt=1,\n",
" nlat=1,\n",
" nlon=1,\n",
" biovert_opt=[0, 1, 2, 3],\n",
")\n",
"ds = run_config_sens(cases)\n",
"ds"
]
},
{
"cell_type": "markdown",
"id": "affd849c-3464-4f27-abc3-1390a4395e31",
"metadata": {},
"source": [
"Currently, method 1 should give the same result as canopy-integrating the laywerwise method 0 results afterwards.\n",
"\n",
"The method 2 result is smaller since it places more weight on the lower--middle regions of the canopy, where light levels are lower. Method 3 is somewhere in between these (consider the shape of the LAD profile vs a Gaussian in height)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "86b197ac-9ecb-40b2-a7e5-1354fcc0a2d3",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"dz = ds.z.diff(\"z\")\n",
"\n",
"e0 = float((ds.emi_isop.isel(case=0) * dz).sum(\"z\"))\n",
"e1 = float(ds.emi_isop.isel(case=1, z=-1))\n",
"e2 = float(ds.emi_isop.isel(case=2, z=-1))\n",
"e3 = float(ds.emi_isop.isel(case=3, z=-1))\n",
"print(e0, e1, e2, e3, sep=\"\\n\")\n",
"\n",
"assert np.isclose(e0, e1)"
]
}
],
"metadata": {
Expand Down
101 changes: 91 additions & 10 deletions src/canopy_bioemi_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,8 @@ module canopy_bioemi_mod

!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, &
TEMP2, LU_OPT, VTYPE, MODRES, CCE, EMI_IND, EMI_OUT)
TEMP2, LU_OPT, VTYPE, MODRES, CCE, VERT, CO2OPT, CO2SET, &
EMI_IND, EMI_OUT)

!-----------------------------------------------------------------------

Expand All @@ -32,7 +33,7 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, &
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
use canopy_const_mod, ONLY: rk,rgasuniv !constants for canopy models
use canopy_utils_mod, ONLY: interp_linear1_internal
use canopy_utils_mod, ONLY: interp_linear1_internal, GET_GAMMA_CO2
use canopy_phot_mod
use canopy_bioparm_mod

Expand All @@ -51,6 +52,9 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, &
INTEGER, INTENT( IN ) :: VTYPE ! Grid cell dominant vegetation type
REAL(RK), INTENT( IN ) :: MODRES ! Canopy model input vertical resolution (m)
REAL(RK), INTENT( IN ) :: CCE ! MEGAN Canopy environment coefficient.
INTEGER, INTENT( IN ) :: VERT ! MEGAN vertical integration option (default = 0/no integration)
INTEGER, INTENT( IN ) :: CO2OPT ! Option for co2 inhibition calculation
REAL(RK), INTENT( IN ) :: CO2SET ! User set atmospheric CO2 conc [ppmv]
INTEGER, INTENT( IN ) :: EMI_IND ! Input biogenic emissions index
REAL(RK), INTENT( OUT ) :: EMI_OUT(:) ! Output canopy layer volume emissions (kg m-3 s-1)

Expand Down Expand Up @@ -92,10 +96,13 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, &
REAL(RK) :: E_OPT(SIZE(ZK)) ! maximum normalized emission capacity
REAL(RK) :: TLEAF_OPT(SIZE(ZK)) ! Tleaf at which E_OPT occurs (K)
REAL(RK) :: FLAI(SIZE(ZK)) ! Fractional LAI in layer
REAL(RK) :: VPGWT(SIZE(ZK)) ! MEGANv3-like in-canopy weighting factor
REAL(RK) :: GAUSS(SIZE(ZK)) ! MEGANv3-like in-canopy gaussian
REAL(RK) :: CT1 ! Activation energy (kJ/mol)
REAL(RK) :: CEO ! Empirical coefficient
REAL(RK) :: EF ! Final Mapped Emission factor (EF) (ug/m2 hr)
integer i
REAL(RK) :: GAMMACO2 ! CO2 inhibition factor (isoprene only)
integer i, LAYERS

! Constant Canopy Parameters
REAL(RK), PARAMETER :: FRAC_PAR = 0.5_rk !Fraction of incoming solar irradiance that is PAR
Expand Down Expand Up @@ -298,15 +305,89 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, CLU, COSZEN, SFCRAD, &

GammaPPFD_AVE = (GammaPPFD_SUN*FSUN) + (GammaPPFD_SHADE*(1.0-FSUN)) ! average = sum sun and shade weighted by sunlit fraction

! Get CO2 inhibition factor for isoprene only

if (EMI_IND .eq. 1) then !Isoprene
GAMMACO2 = GET_GAMMA_CO2(CO2OPT,CO2SET)
else
GAMMACO2 = 1.0_rk
end if

! Calculate emissions profile in the canopy
EMI_OUT = 0.0_rk ! set initial emissions profile to zero
do i=1, SIZE(ZK)
if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then ! above ground level and at/below canopy top
FLAI(i) = ((FCLAI(i+1) - FCLAI(i)) * LAI)/MODRES !fractional LAI in each layer converted to LAD (m2 m-3)
EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * CCE ! (ug m-3 hr-1)
EMI_OUT(i) = EMI_OUT(i) * 2.7777777777778E-13_rk !TBD: convert emissions output to (kg m-3 s-1)
end if
end do
FLAI = 0.0_rk ! set initial fractional FLAI (LAD) profile to zero

if (VERT .eq. 0) then !Full 3D leaf-level biogenic emissions (no averaging, summing, or integration)
do i=1, SIZE(ZK)
if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then ! above ground level and at/below canopy top
FLAI(i) = ((FCLAI(i+1) - FCLAI(i)) * LAI)/MODRES !fractional LAI in each layer converted to LAD (m2 m-3)
EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * GAMMACO2 * CCE ! (ug m-3 hr-1)
EMI_OUT(i) = EMI_OUT(i) * 2.7777777777778E-13_rk !convert emissions output to (kg m-3 s-1)
end if
end do
else if (VERT .eq. 1) then !"MEGANv3-like": Use weighting factors normalized to plant distribution shape (FCLAI)
!across canopy layers
LAYERS = floor(FCH/MODRES) + 1
do i=1, SIZE(ZK)
if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then
FLAI(i) = ((FCLAI(i+1) - FCLAI(i)) * LAI)/MODRES !fractional LAI in each layer converted to LAD (m2 m-3)
end if
end do
do i=1, SIZE(ZK)
if (ZK(i) .gt. 0.0 .and. ZK(i) .le. FCH) then
VPGWT(i) = (FLAI(i))/sum(FLAI(1:LAYERS))
end if
end do
EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * &
VPGWT(1:LAYERS)) * GAMMACO2 * CCE !put into top model layer (ug m-2 hr-1)
EMI_OUT = EMI_OUT * 2.7777777777778E-13_rk !convert emissions output to (kg m-2 s-1)
else if (VERT .eq. 2) then !"MEGANv3-like": Add weighted sum of activity coefficients using normal distribution
!across canopy layers using 5 layer numbers directly from MEGANv3
!--warning: weights are not consistent with FCLAI distribution
!used for biomass distribution used for sunlit/shaded in Gamma TLEAF and GammaPPFD.
LAYERS = floor(FCH/MODRES) + 1
do i=1, SIZE(ZK)
if (ZK(i) .gt. FCH) then
GAUSS(i) = 0.0
else if (ZK(i) .le. FCH .and. ZK(i) .gt. FCH*(4.0_rk/5.0_rk)) then !Level 1 - 2
GAUSS(i) = interp_linear1_internal((/ FCH*(4.0_rk/5.0_rk),FCH /), &
(/ 0.118464_rk,0.0_rk /),ZK(i))
else if (ZK(i) .le. FCH*(4.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(3.0_rk/5.0_rk)) then !Level 2 - 3
GAUSS(i) = interp_linear1_internal((/ FCH*(3.0_rk/5.0_rk),FCH*(4.0_rk/5.0_rk) /), &
(/ 0.239314_rk,0.118464_rk /),ZK(i))
else if (ZK(i) .le. FCH*(3.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(2.0_rk/5.0_rk)) then !Level 3 - 4
GAUSS(i) = interp_linear1_internal((/ FCH*(2.0_rk/5.0_rk),FCH*(3.0_rk/5.0_rk) /), &
(/ 0.284444_rk,0.239314_rk /),ZK(i))
else if (ZK(i) .le. FCH*(2.0_rk/5.0_rk) .and. ZK(i) .gt. FCH*(1.0_rk/5.0_rk) ) then !Level 4 - Bottom
GAUSS(i) = interp_linear1_internal((/ FCH*(1.0_rk/5.0_rk),FCH*(2.0_rk/5.0_rk) /), &
(/ 0.239314_rk,0.284444_rk /),ZK(i))
else if (ZK(i) .le. FCH*(1.0_rk/5.0_rk) ) then !Level 4 - Bottom
GAUSS(i) = interp_linear1_internal((/ ZK(1),FCH*(1.0_rk/5.0_rk) /), &
(/ 0.118464_rk,0.239314_rk /),ZK(i))
end if
end do

do i=1, SIZE(ZK)
VPGWT(i) = GAUSS(i)/sum(GAUSS(1:LAYERS))
end do
EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * &
VPGWT(1:LAYERS)) * GAMMACO2 * CCE !put into top model layer (ug m-2 hr-1)
EMI_OUT = EMI_OUT * 2.7777777777778E-13_rk !convert emissions output to (kg m-2 s-1)
else if (VERT .eq. 3) then !"MEGANv3-like": Add weighted sum of activity coefficients equally
!across canopy layers
!--warning: weights are not consistent with FCLAI distribution
!used for biomass distribution used for sunlit/shaded in Gamma TLEAF and GammaPPFD.
LAYERS = floor(FCH/MODRES) + 1
do i=1, SIZE(ZK)
VPGWT(i) = 1.0_rk/LAYERS
end do
EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * &
VPGWT(1:LAYERS)) * GAMMACO2 * CCE !put into top model layer (ug m-2 hr-1)
EMI_OUT = EMI_OUT * 2.7777777777778E-13_rk !convert emissions output to (kg m-2 s-1)
else
write(*,*) 'Wrong BIOVERT_OPT choice of ', VERT, ' in namelist...exiting'
call exit(2)
end if

END SUBROUTINE CANOPY_BIO

Expand Down
Loading

0 comments on commit 1a5dabc

Please sign in to comment.