diff --git a/README.md b/README.md index 409f0e3e..f3ba4ed7 100644 --- a/README.md +++ b/README.md @@ -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) | diff --git a/input/namelist.canopy b/input/namelist.canopy index a07e3417..04781010 100755 --- a/input/namelist.canopy +++ b/input/namelist.canopy @@ -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 / diff --git a/python/examples.ipynb b/python/examples.ipynb index d9a75e91..62310959 100644 --- a/python/examples.ipynb +++ b/python/examples.ipynb @@ -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": { diff --git a/src/canopy_bioemi_mod.F90 b/src/canopy_bioemi_mod.F90 index 85da7d7c..ba3d87d2 100644 --- a/src/canopy_bioemi_mod.F90 +++ b/src/canopy_bioemi_mod.F90 @@ -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) !----------------------------------------------------------------------- @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/src/canopy_calcs.F90 b/src/canopy_calcs.F90 index 3f47929b..e4aef172 100644 --- a/src/canopy_calcs.F90 +++ b/src/canopy_calcs.F90 @@ -84,6 +84,38 @@ SUBROUTINE canopy_calcs if (lu_opt .eq. 0 .or. lu_opt .eq. 1 ) then !VIIRS or MODIS if (vtyperef .gt. 0 .and. vtyperef .le. 10 .or. vtyperef .eq. 12) then !VIIRS or MODIS types +! ... check for ssg_opt from user namelist + if (vtyperef .ge. 6 .and. vtyperef .le. 10) then !VIIRS/MODIS shrubs/savannas/grasses (SSG) type + if (ssg_opt .eq. 0) then !use GEDI inputs for SSG heights (not likely captured...) + hcmref = hcmref + else if (ssg_opt .eq. 1) then !user set constant shrubs/savannas/grasslands height + hcmref = ssg_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong SSG_OPT choice of ', ssg_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + +! ... check for crop_opt from user namelist + if (vtyperef .eq. 12) then !VIIRS/MODIS crop type + if (crop_opt .eq. 0) then !use GEDI inputs for crop height (not likely captured...) + hcmref = hcmref + else if (crop_opt .eq. 1) then !user set constant crop height + hcmref = crop_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong CROP_OPT choice of ', crop_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + ! ... check for contiguous canopy conditions at each model grid cell if (hcmref .gt. fch_thresh .and. ffracref .gt. frt_thresh & .and. lairef .gt. lai_thresh) then @@ -153,79 +185,98 @@ SUBROUTINE canopy_calcs !ISOP call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 1, emi_isop_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 1, emi_isop_3d(i,j,:)) !MYRC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 2, emi_myrc_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 2, emi_myrc_3d(i,j,:)) !SABI call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 3, emi_sabi_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 3, emi_sabi_3d(i,j,:)) !LIMO call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 4, emi_limo_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 4, emi_limo_3d(i,j,:)) !CARE call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 5, emi_care_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 5, emi_care_3d(i,j,:)) !OCIM call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 6, emi_ocim_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 6, emi_ocim_3d(i,j,:)) !BPIN call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 7, emi_bpin_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 7, emi_bpin_3d(i,j,:)) !APIN call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 8, emi_apin_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 8, emi_apin_3d(i,j,:)) !MONO call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 9, emi_mono_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 9, emi_mono_3d(i,j,:)) !FARN call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 10, emi_farn_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 10, emi_farn_3d(i,j,:)) !CARY call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 11, emi_cary_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 11, emi_cary_3d(i,j,:)) !SESQ call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 12, emi_sesq_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 12, emi_sesq_3d(i,j,:)) !MBOL call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 13, emi_mbol_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 13, emi_mbol_3d(i,j,:)) !METH call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 14, emi_meth_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 14, emi_meth_3d(i,j,:)) !ACET call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 15, emi_acet_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 15, emi_acet_3d(i,j,:)) !CO call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 16, emi_co_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 16, emi_co_3d(i,j,:)) !BIDI VOC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 17, emi_bvoc_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 17, emi_bvoc_3d(i,j,:)) !Stress VOC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 18, emi_svoc_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 18, emi_svoc_3d(i,j,:)) !Other VOC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 19, emi_ovoc_3d(i,j,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 19, emi_ovoc_3d(i,j,:)) end if end if @@ -296,6 +347,38 @@ SUBROUTINE canopy_calcs if (lu_opt .eq. 0 .or. lu_opt .eq. 1 ) then !VIIRS or MODIS if (vtyperef .gt. 0 .and. vtyperef .le. 10 .or. vtyperef .eq. 12) then !VIIRS or MODIS types +! ... check for ssg_opt from user namelist + if (vtyperef .ge. 6 .and. vtyperef .le. 10) then !VIIRS/MODIS shrubs/savannas/grasses (SSG) type + if (ssg_opt .eq. 0) then !use GEDI inputs for SSG heights (not likely captured...) + hcmref = hcmref + else if (ssg_opt .eq. 1) then !user set constant shrubs/savannas/grasslands height + hcmref = ssg_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong SSG_OPT choice of ', ssg_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + +! ... check for crop_opt from user namelist + if (vtyperef .eq. 12) then !VIIRS/MODIS crop type + if (crop_opt .eq. 0) then !use GEDI inputs for crop height + hcmref = hcmref + else if (crop_opt .eq. 1) then !user set constant crop height + hcmref = crop_set + !recalculate + zhc = zk/hcmref + cansublays = floor(hcmref/modres) + else + write(*,*) 'Wrong CROP_OPT choice of ', crop_opt, & + ' in namelist...exiting' + call exit(2) + end if + end if + ! ... check for contiguous canopy conditions at each model grid cell if (hcmref .gt. fch_thresh .and. ffracref .gt. frt_thresh & .and. lairef .gt. lai_thresh) then @@ -365,79 +448,98 @@ SUBROUTINE canopy_calcs !ISOP call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 1, emi_isop(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 1, emi_isop(loc,:)) !MYRC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 2, emi_myrc(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 2, emi_myrc(loc,:)) !SABI call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 3, emi_sabi(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 3, emi_sabi(loc,:)) !LIMO call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 4, emi_limo(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 4, emi_limo(loc,:)) !CARE call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 5, emi_care(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 5, emi_care(loc,:)) !OCIM call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 6, emi_ocim(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 6, emi_ocim(loc,:)) !BPIN call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 7, emi_bpin(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 7, emi_bpin(loc,:)) !APIN call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 8, emi_apin(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 8, emi_apin(loc,:)) !MONO call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 9, emi_mono(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 9, emi_mono(loc,:)) !FARN call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 10, emi_farn(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 10, emi_farn(loc,:)) !CARY call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 11, emi_cary(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 11, emi_cary(loc,:)) !SESQ call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 12, emi_sesq(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 12, emi_sesq(loc,:)) !MBOL call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 13, emi_mbol(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 13, emi_mbol(loc,:)) !METH call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 14, emi_meth(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 14, emi_meth(loc,:)) !ACET call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 15, emi_acet(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 15, emi_acet(loc,:)) !CO call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 16, emi_co(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 16, emi_co(loc,:)) !BIDI VOC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 17, emi_bvoc(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 17, emi_bvoc(loc,:)) !Stress VOC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 18, emi_svoc(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 18, emi_svoc(loc,:)) !Other VOC call canopy_bio(zk, fafraczInt, hcmref, & lairef, cluref, cszref, dswrfref, tmp2mref, & - lu_opt, vtyperef, modres, bio_cce, 19, emi_ovoc(loc,:)) + lu_opt, vtyperef, modres, bio_cce, biovert_opt, co2_opt, co2_set, & + 19, emi_ovoc(loc,:)) end if end if diff --git a/src/canopy_canopts_mod.F90 b/src/canopy_canopts_mod.F90 index 2fd44171..b3f1c29c 100644 --- a/src/canopy_canopts_mod.F90 +++ b/src/canopy_canopts_mod.F90 @@ -33,5 +33,12 @@ MODULE canopy_canopts_mod real(rk) :: z0ghc !ratio of ground roughness length to canopy top height real(rk) :: lambdars !Value representing influence of roughness sublayer (nondimensional) real(rk) :: bio_cce !MEGAN biogenic emission canopy environment coefficient. + integer :: biovert_opt !MEGAN vertical integration of emissions option (default = 0, off) + integer :: ssg_opt !Set default integer for shrubs/savanna/grassland vegtype option from GEDI or user (default = 0) + real(rk) :: ssg_set !Set default value for shrubs/savanna/grassland vegtype heights used in model (m) (Default = 1 m) + integer :: crop_opt !Set default integer for crop vegtype option from GEDI or user (default = 0) + real(rk) :: crop_set !Set default value for crop vegtype heights used in model (m) (Default = 3 m) + integer :: co2_opt !Set default integer for co2 inhibition option for biogenic isoprene emissions (default=0; Possell & Hewitt (2011)) + real(rk) :: co2_set !Set default value for atmospheric co2 concentration for co2_opt (m) (Default = 400.0 ppmv) END MODULE canopy_canopts_mod diff --git a/src/canopy_readnml.F90 b/src/canopy_readnml.F90 index a9eb4f4a..52e63989 100644 --- a/src/canopy_readnml.F90 +++ b/src/canopy_readnml.F90 @@ -23,7 +23,8 @@ SUBROUTINE canopy_readnml NAMELIST /userdefs/ infmt_opt, nlat, nlon, modlays, modres, href_opt, href_set, & z0ghc, lambdars, flameh_opt, flameh_set, frp_fac, ifcanwind, ifcanwaf, & ifcaneddy, ifcanphot, ifcanbio, pai_opt, pai_set, lu_opt, z0_opt, dx_opt, & - dx_set, lai_thresh, frt_thresh, fch_thresh, rsl_opt, bio_cce + dx_set, lai_thresh, frt_thresh, fch_thresh, rsl_opt, bio_cce, biovert_opt, & + ssg_opt, ssg_set, crop_opt, crop_set, co2_opt, co2_set !------------------------------------------------------------------------------- ! Error, warning, and informational messages. @@ -201,6 +202,41 @@ SUBROUTINE canopy_readnml bio_cce = 0.21_rk !------------------------------------------------------------------------------- +!------------------------------------------------------------------------------- +! Set default integer value for MEGAN vertical integration of emissions (0, off full leaf-level emissions) + biovert_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default integer for shrubs/savanaa/grasslands vegtype option from GEDI or user (default = 0) + ssg_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default value for shrubs/savanaa/grasslands vegtype heights used in model (m) (Default = 1 m) + ssg_set = 1.0_rk +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default integer for crop vegtype option from GEDI or user (default = 0) + crop_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default value for crop vegtype heights used in model (m) (Default = 3 m) + crop_set = 3.0_rk +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default integer for co2 inhibition option for biogenic isoprene emissions (default = 0; Possell & Hewitt (2011)) + co2_opt = 0 +!------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +! Set default value for atmospheric co2 concentration for co2_opt (m) (Default = 400.0 ppmv) + co2_set = 400.0_rk +!------------------------------------------------------------------------------- + !------------------------------------------------------------------------------- ! Read namelist to get user definitions. Rewind namelist file after each ! read in case namelists are not in the correct order in the namelist. diff --git a/src/canopy_utils_mod.F90 b/src/canopy_utils_mod.F90 index 312f07a3..ae1117bb 100644 --- a/src/canopy_utils_mod.F90 +++ b/src/canopy_utils_mod.F90 @@ -6,7 +6,7 @@ module canopy_utils_mod private public IntegrateTrapezoid,interp_linear1_internal,CalcPAI, & - CalcDX,CalcFlameH + CalcDX,CalcFlameH,GET_GAMMA_CO2 contains @@ -125,4 +125,130 @@ function CalcFlameH(frp, dx) end function + real(rk) function GET_GAMMA_CO2(co2_opt, co2_set) result( GAMMA_CO2 ) + ! !IROUTINE: get_gamma_co2 + ! + ! !DESCRIPTION: Function GET\_GAMMA\_CO2 computes the CO2 activity factor + ! associated with CO2 inhibition of isoprene emission. Called from + ! GET\_MEGAN\_EMISSIONS only. + !\\ + !\\ + ! !INTERFACE: + + ! !INPUT PARAMETERS: + integer, INTENT(IN) :: co2_opt ! Option for co2 inhibition calculation + ! 0=Possell & Hewitt (2011); + ! 1=Wilkinson et al. (2009) + ! >1=off + real(rk), INTENT(IN) :: co2_set ! User set atmospheric CO2 conc [ppmv] + ! + ! !RETURN VALUE: +! REAL(rk) :: GAMMA_CO2 ! CO2 activity factor [unitless] + ! + ! !LOCAL VARIABLES: + REAL(rk) :: CO2i ! Intercellular CO2 conc [ppmv] + REAL(rk) :: ISMAXi ! Asymptote for intercellular CO2 + REAL(rk) :: HEXPi ! Exponent for intercellular CO2 + REAL(rk) :: CSTARi ! Scaling coef for intercellular CO2 + REAL(rk) :: ISMAXa ! Asymptote for atmospheric CO2 + REAL(rk) :: HEXPa ! Exponent for atmospheric CO2 + REAL(rk) :: CSTARa ! Scaling coef for atmospheric CO2 + ! + ! !REMARKS: + ! References: + ! ============================================================================ + ! (1 ) Heald, C. L., Wilkinson, M. J., Monson, R. K., Alo, C. A., + ! Wang, G. L., and Guenther, A.: Response of isoprene emission + ! to ambient co(2) changes and implications for global budgets, + ! Global Change Biology, 15, 1127-1140, 2009. + ! (2 ) Wilkinson, M. J., Monson, R. K., Trahan, N., Lee, S., Brown, E., + ! Jackson, R. B., Polley, H. W., Fay, P. A., and Fall, R.: Leaf + ! isoprene emission rate as a function of atmospheric CO2 + ! concentration, Global Change Biology, 15, 1189-1200, 2009. + ! (3 ) Possell, M., and Hewitt, C. N.: Isoprene emissions from plants + ! are mediated by atmospheric co2 concentrations, Global Change + ! Biology, 17, 1595-1610, 2011. + + ! !REVISION HISTORY: + ! (1 ) Implemented in the standard code by A. Tai (Jun 2012). + ! See https://github.com/geoschem/hemco for complete history + !EOP + !------------------------------------------------------------------------------ + !BOC + + !----------------------- + ! Compute GAMMA_CO2 + !----------------------- + + !---------------------------------------------------------- + ! Choose between two alternative CO2 inhibition schemes + !---------------------------------------------------------- + + IF ( co2_opt .eq. 0 ) THEN + + ! Use empirical relationship of Possell & Hewitt (2011): + ! Empirical relationship of Possell & Hewitt (2011) based on nine + ! experimental studies including Wilkinson et al. (2009). + + GAMMA_CO2 = 8.9406_rk / ( 1.0_rk + 8.9406_rk * 0.0024_rk * co2_set ) + + + ELSEIF ( co2_opt .eq. 1 ) THEN + + ! Use parameterization of Wilkinson et al. (2009): + ! Semi-process-based parameterization of Wilkinson et al. (2009), + ! taking into account of sensitivity to intercellular CO2 + ! fluctuation, which is here set as a constant fraction of + ! atmospheric CO2. This is especially recommended for sub-ambient + ! CO2 concentrations:: + + + ! Parameters for intercellular CO2 using linear interpolation: + IF ( co2_set <= 600.0_rk ) THEN + ISMAXi = 1.036_rk - (1.036_rk - 1.072_rk) / & + (600.0_rk - 400.0_rk) * (600.0_rk - co2_set) + HEXPi = 2.0125_rk - (2.0125_rk - 1.7000_rk) / & + (600.0_rk - 400.0_rk) * (600.0_rk - co2_set) + CSTARi = 1150.0_rk - (1150.0_rk - 1218.0_rk) / & + (600.0_rk - 400.0_rk) * (600.0_rk - co2_set) + ELSEIF ( co2_set > 600.0_rk .AND. co2_set < 800.0_rk ) THEN + ISMAXi = 1.046_rk - (1.046_rk - 1.036_rk) / & + (800.0_rk - 600.0_rk) * (800.0_rk - co2_set) + HEXPi = 1.5380_rk - (1.5380_rk - 2.0125_rk) / & + (800.0_rk - 600.0_rk) * (800.0_rk - co2_set) + CSTARi = 2025.0_rk - (2025.0_rk - 1150.0_rk) / & + (800.0_rk - 600.0_rk) * (800.0_rk - co2_set) + ELSE + ISMAXi = 1.014_rk - (1.014_rk - 1.046_rk) / & + (1200.0_rk - 800.0_rk) * (1200.0_rk - co2_set) + HEXPi = 2.8610_rk - (2.8610_rk - 1.5380_rk) / & + (1200.0_rk - 800.0_rk) * (1200.0_rk - co2_set) + CSTARi = 1525.0_rk - (1525.0_rk - 2025.0_rk) / & + (1200.0_rk - 800.0_rk) * (1200.0_rk - co2_set) + ENDIF + + ! Parameters for atmospheric CO2: + ISMAXa = 1.344_rk + HEXPa = 1.4614_rk + CSTARa = 585.0_rk + + ! For now, set CO2_Ci = 0.7d0 * CO2_Ca as recommended by Heald + ! et al. (2009): + CO2i = 0.7_rk * co2_set + + ! Compute GAMMA_CO2: + GAMMA_CO2 = ( ISMAXi - ISMAXi * CO2i**HEXPi / & + ( CSTARi**HEXPi + CO2i**HEXPi ) ) & + * ( ISMAXa - ISMAXa * ( 0.7_rk * co2_set )**HEXPa / & + ( CSTARa**HEXPa + ( 0.7_rk * co2_set )**HEXPa ) ) + + ELSE + + ! No CO2 inhibition scheme is used; GAMMA_CO2 set to unity: + GAMMA_CO2 = 1.0_rk + + ENDIF + + end function GET_GAMMA_CO2 + end module canopy_utils_mod