Adding Gamma AQ, HT, LT, and HW stress factors to biogenics

README.md

@@ -85,6 +85,13 @@ Current Canopy-App components:
 
 ## Outputs
 
+Namelist Option : `file_out`  Prefix string (e.g., `'test'`) used to name output file (Output is 1D txt when using input 1D data (i.e., `infmt_opt=1`), or is 2D NetCDF output when 2D NetCDF input is used (i.e., `infmt_opt=0`)).
+
+Current 3D fields include canopy winds (`canwind`), canopy vertical/eddy diffusivity values `kz`), biogenic emissions (see Table 1 below),
+canopy photolysis attenuation correction factors (`rjcf`), and derived Leaf Area Density (`lad`) from the foliage shape function.  
+
+Current 2D fields includes the Wind Adjustment Factor (`waf`), flame heights (`flameh`), and canopy heights (`canheight`). Current 1D fields include the canopy model interface levels (`z`).
+
 **Note for Biogenic emissions:** When `ifcanbio=.TRUE.`, output will include 3D canopy resolved biogenic emissions for the following species (based on Guenther et al., 2012), which have been mapped from Guenther et al. PFTs to input LU_OPT.
 
 ### Table 1. Canopy-App Biogenic Emissions Output Variables
@@ -111,11 +118,6 @@ Current Canopy-App components:
 | `emi_svoc`    | Stress VOC emissions (15 compounds, Table 1 Guenther et al. (2012) | 18        |
 | `emi_ovoc`    | Other VOC emissions (49 compounds, Table 1 Guenther et al. (2012) | 19        |
 
-**Current Canopy-App Output:** As discussed above, the current Canopy-App optional outputs includes 3D canopy winds (`canwind`), canopy vertical/eddy diffusivity values `kz`), biogenic emissions (see Table 1), and
-canopy photolysis attenuation correction factors (`rjcf`), and derived Leaf Area Density (`lad`) from the foliage shape function.  Current 2D fields includes the Wind Adjustment Factor (`waf`), flame heights (`flameh`), and canopy heights (`canheight`). Current 1D fields include the canopy model interface levels (`z`).
-
-Namelist Option : `file_out`  Prefix string (e.g., `'test'`) used to name output file (Output is 1D txt when using input 1D data (i.e., `infmt_opt=1`), or is 2D NetCDF output when 2D NetCDF input is used (i.e., `infmt_opt=0`)).
-
 ## Inputs and Settings
 
 **Current Canopy-App Input:** Typical 1D or 2D (time=1,lat,lon) gridded atmospheric model input variables in 1st layer above canopy.  Some 3D inputs are supported (see `var3d_opt` in [Table 3](#table-3-current-user-namelist-options) and associated options).
@@ -170,7 +172,7 @@ The Canopy-App input data in [Table 2](#table-2-canopy-app-required-input-variab
 | `csz`                            | Cosine of the solar zenith angle (dimensionless) | [Based on Python Pysolar](https://pysolar.readthedocs.io/en/latest/) |
 | `mol`                            | Monin-Obukhov Length (m)                    | Externally calculated using GFS `tmp2m`, `fricv`, and `shtfl`.  ([Essa, 1999](https://inis.iaea.org/collection/NCLCollectionStore/_Public/37/118/37118528.pdf)) |
 | `href`                           | Reference height above canopy (m) - 10 m    | Assumed constant (i.e., 10 m).  Can be taken from NL. |
-| `ozone_w126`                     | Ozone W126 index (ppm-hours)                | Climatological value during 04/2021-04/2024 based on GFSv16 lowest model layer ozone mixing ratios. |
+| `ozone_w126`                     | Ozone W126 index (ppm-hours)                | A three year climatological calculation between 04/2021-04/2024, based on GFSv16 lowest model layer ozone mixing ratios.  The W126 calculation is based on the [EPA definition](https://www.epa.gov/sites/default/files/2015-09/documents/w126_steps_to_calculate_revised_feb19.pdf). |
 
 **More Information on Data Sources from [Table 2](#table-2-canopy-app-required-input-variables):**
 
@@ -282,6 +284,12 @@ Otherwise, please contact [email protected] for other GFSv16 data peri
 | `hist_opt`      | user-set option to use historically averaged short-term (24-hr) and long-term (240-hr) rolling averages for leaf temperature and PAR for biogenic emissions  (default is off i.e., `hist_opt=0`)  Note: If simulation is </= 24 hours, instantaneous values for leaf temperature and PAR will be used even if historical averaging is turned on (i.e., `hist_opt=1`).  Recommend turning on `hist_opt=1` and running at least for 25 hours to create a model spin-up, and use subsequent simulation hours for analysis. Overall a 10-day (240 hr) spinup is optimal for best analysis of biogenic emissions. |
 | `soim_opt`   | user-set options for applying soil moisture response to biogenic VOC emissions based on [Guenther et al. 2006](https://doi.org/10.5194/acp-6-3181-2006), which depends on input soil moisture at depth and the wilting point.  This includes additional PFT dependent approach for cumulative root fraction within each soil layer from [Zeng (2001)](https://doi.org/10.1175/1525-7541(2001)002<0525:GVRDFL>2.0.CO;2). (default is off i.e., `soim_opt=1`, the corresponding $\gamma$ is set to 1). If turned on (`soim_opt=0`), which is recommended, soim $\gamma$ is calculated and the prescribed 4-layer soil depths (`soild[1-4]` below) are used.  Four layers are assumed, and are based on GFS Noah and Noah-MP LSM. |
 | `soild[1-4]` | user-set real values of four level soil depths at centerpoint (cm).  Four layers are based on the GFS Noah and Noah-MP LSM, default values are `soild1=5.0`, `soild2=25.0`, `soild3=70.0`, and  `soild4=150.0`. |
+| `aq_opt`       | user-set options for applying an air quality stress factor for biogenic emissions using calculated, spatially-dependent and global GFS-based ozone W126 values (= `0`) or a constant user-set W126 value (= `1`). To turn off aq stress factor set `aq_opt=2` (set as default, Off).  Note:  The aq_opt should only be turned on during simulations across respective ozone season for specific region (e.g., April-October in the U.S.)|
+| `w126_set`       | user-set real value of constant ozone W126 values (ppm-hours) (only used if `aq_opt=1`) |
+| `ht_opt`       | user-set options for applying a daily high temperature stress factor for biogenic emissions using daily maximum 2-meter input temperature (= `0`).  This is based on MEGAN3 and it is recommended that this option is only used when turning the historical option on (i.e., `hist_opt=1`) and running longer than 1-day simulations to obtain the daily max.   To turn off ht stress factor set `ht_opt=1` (set as default, Off)  |
+| `lt_opt`       | user-set options for applying a daily low temperature stress factor for biogenic emissions using daily minimum 2-meter input temperature (= `0`).  This is based on MEGAN3 and it is recommended that this option is only used when turning the historical option on (i.e., `hist_opt=1`) and running longer than 1-day simulations to obtain the daily min.   To turn off lt stress factor set `lt_opt=1` (set as default, Off) |
+| `hw_opt`       | user-set options for applying a daily high wind speed stress factor for biogenic emissions using daily maximum 10-meter input wind speed (= `0`).  This is based on MEGAN3 and it is recommended that this option is only used when turning the historical option on (i.e., `hist_opt=1`) and running longer than 1-day simulations to obtain the daily max.   To turn off hw stress factor set `hw_opt=1` (set as default, Off) |
+
 
 **\*\*** If `modres` >> `flameh` then some error in WAF calculation will be incurred.  Suggestion is to use relative fine `modres` (at least <= 0.5 m) compared to average flame heights (e.g., ~ 1.0 m) if WAF is required.
 

input/namelist.canopy

@@ -83,5 +83,10 @@
   soild2      =  25.0
   soild3      =  70.0
   soild4      =  150.0
+  aq_opt      =  2
+  w126_set    =  30.0
   hist_opt    =  1
+  ht_opt      =  1
+  lt_opt      =  1
+  hw_opt      =  1
 /

src/canopy_alloc.F90

@@ -67,6 +67,12 @@ SUBROUTINE canopy_alloc
     if(.not.allocated(ppfd_sun240))        allocate(ppfd_sun240(nlat*nlon,modlays))
     if(.not.allocated(ppfd_shade240))      allocate(ppfd_shade240(nlat*nlon,modlays))
 
+    if(.not.allocated(tmp2mref_tmp))       allocate(tmp2mref_tmp(ntime,nlat*nlon))
+    if(.not.allocated(ubzref_tmp))         allocate(ubzref_tmp(ntime,nlat*nlon))
+    if(.not.allocated(daily_maxt2m))       allocate(daily_maxt2m(nlat*nlon))
+    if(.not.allocated(daily_mint2m))       allocate(daily_mint2m(nlat*nlon))
+    if(.not.allocated(daily_maxws10m))     allocate(daily_maxws10m(nlat*nlon))
+
     if(.not.allocated(tleaf_sun24_tmp_3d))    allocate(tleaf_sun24_tmp_3d(ntime,nlon,nlat,modlays))
     if(.not.allocated(tleaf_shade24_tmp_3d))  allocate(tleaf_shade24_tmp_3d(ntime,nlon,nlat,modlays))
     if(.not.allocated(tleaf_ave24_tmp_3d))    allocate(tleaf_ave24_tmp_3d(ntime,nlon,nlat,modlays))
@@ -88,6 +94,12 @@ SUBROUTINE canopy_alloc
     if(.not.allocated(ppfd_sun240_3d))        allocate(ppfd_sun240_3d(nlon,nlat,modlays))
     if(.not.allocated(ppfd_shade240_3d))      allocate(ppfd_shade240_3d(nlon,nlat,modlays))
 
+    if(.not.allocated(tmp2mref_tmp_3d))       allocate(tmp2mref_tmp_3d(ntime,nlon,nlat))
+    if(.not.allocated(ubzref_tmp_3d))         allocate(ubzref_tmp_3d(ntime,nlon,nlat))
+    if(.not.allocated(daily_maxt2m_2d))       allocate(daily_maxt2m_2d(nlon,nlat))
+    if(.not.allocated(daily_mint2m_2d))       allocate(daily_mint2m_2d(nlon,nlat))
+    if(.not.allocated(daily_maxws10m_2d))     allocate(daily_maxws10m_2d(nlon,nlat))
+
 !-------------------------------------------------------------------------------
 ! Allocate arrays for Canopy Wind Outputs
 !-------------------------------------------------------------------------------

src/canopy_bioemi_mod.F90

@@ -14,7 +14,9 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
         LEAFAGEOPT, PASTLAI, CURRENTLAI, TSTEPLAI, &
         LOSSOPT, LOSSSET, LOSSIND, LIFETIME, USTAR, &
         SOIMOPT, SOIM1, SOIM2, SOIM3, SOIM4, SOID1, SOID2, SOID3, &
-        SOID4, WILT, &
+        SOID4, WILT, AQOPT, W126_SET, W126_REF,  &
+        HTOPT, LTOPT, HWOPT, DAILY_MAXT2, DAILY_MINT2, &
+        DAILY_MAXWS10, &
         MODLAYS, EMI_IND, EMI_OUT)
 
 !-----------------------------------------------------------------------
@@ -47,7 +49,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
         use canopy_const_mod,  ONLY: rk,rgasuniv   !constants for canopy models
         use canopy_utils_mod,  ONLY: interp_linear1_internal, &
             GET_GAMMA_CO2,GET_GAMMA_LEAFAGE, &
-            GET_GAMMA_SOIM, GET_CANLOSS_BIO
+            GET_GAMMA_SOIM, GET_GAMMA_AQ, GET_GAMMA_HT, &
+            GET_GAMMA_LT, GET_GAMMA_HW, GET_CANLOSS_BIO
         use canopy_bioparm_mod
         use canopy_tleaf_mod
 
@@ -88,6 +91,15 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
         REAL(RK),    INTENT( IN )       :: SOID3           ! Soil depth layer 3 [cm]
         REAL(RK),    INTENT( IN )       :: SOID4           ! Soil depth layer 4 [cm]
         REAL(RK),    INTENT( IN )       :: WILT            ! Wilting point [proportion]
+        INTEGER,     INTENT( IN )       :: AQOPT           ! Option for aq stress calculation
+        REAL(RK),    INTENT( IN )       :: W126_SET        ! User set ozone W126 [ppm-hours]
+        REAL(RK),    INTENT( IN )       :: W126_REF        ! GFS calculated, ozone W126 [ppm-hours]
+        INTEGER,     INTENT( IN )       :: HTOPT           ! Option for high temperature stress calculation
+        INTEGER,     INTENT( IN )       :: LTOPT           ! Option for low temperature stress calculation
+        INTEGER,     INTENT( IN )       :: HWOPT           ! Option for high wind speed stress calculation
+        REAL(RK),    INTENT( IN )       :: DAILY_MAXT2     ! Daily maximum model input 2-m temperature (K)
+        REAL(RK),    INTENT( IN )       :: DAILY_MINT2     ! Daily minimum model input 2-m temperature (K)
+        REAL(RK),    INTENT( IN )       :: DAILY_MAXWS10   ! Daily maximum model input 10-m wind speed temperature (m/s)
 
         INTEGER,    INTENT( IN )        :: LEAFAGEOPT     ! leafage_opt (0= ON, 1= off i.e. GAMMALEAFAGE =1, in canopy_readnml.F90)
         REAL(RK),    INTENT( IN )       :: PASTLAI        ! Past LAI [cm2/cm2]
@@ -137,15 +149,43 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
         REAL(RK) :: AMAT
         REAL(RK) :: AOLD
 
+        !Coeff.'s and threshold/delta threshold for air quality stress factors from canopy_biop
+        REAL(RK) :: CAQ
+        REAL(RK) :: TAQ  ![ppm-hours]
+        REAL(RK) :: DTAQ ![ppm-hours]
+
+        !Coeff.'s and threshold/delta threshold for high temperature stress factors from canopy_biop
+        REAL(RK) :: CHT
+        REAL(RK) :: THT  ![K]
+        REAL(RK) :: DTHT ![K]
+
+        !Coeff.'s and threshold/delta threshold for low temperature stress factors from canopy_biop
+        REAL(RK) :: CLT
+        REAL(RK) :: TLT  ![K]
+        REAL(RK) :: DTLT ![K]
+
+        !Coeff.'s and threshold/delta threshold for high wind stress factors from canopy_biop
+        REAL(RK) :: CHW
+        REAL(RK) :: THW  ![m/s]
+        REAL(RK) :: DTHW ![m/s]
+
+
         ! Coefficients A and B used for PFT dependent cumulative root depth fraction
         REAL(RK) :: ROOTA ! [m-1]
         REAL(RK) :: ROOTB ! [m-1]
         REAL(RK) :: GAMMASOIM ! Soil moisture factor
 
+        REAL(RK) :: GAMMAAQ                        !Air quality stress factor
+        REAL(RK) :: GAMMAHT                        !High temperature stress factor
+        REAL(RK) :: GAMMALT                        !Low temperature stress factor
+        REAL(RK) :: GAMMAHW                        !High wind speed stress factor
+
         REAL(RK) :: GAMMACO2                       ! CO2 inhibition factor (isoprene only)
 
         REAL(RK) :: GAMMALEAFAGE !(SIZE(ZK))                 ! LEAF AGE factor
+
         REAL(RK) :: CANLOSS_FAC                    !Canopy loss factor for summing option
+
         integer i, LAYERS
 
 ! Constant Canopy Parameters
@@ -157,8 +197,12 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
         TLEAF_OPT = 313.0_rk + (0.6_rk * (TLEAF240_AVE-297.0_rk)) !Guenther et al. (2012)
 
 ! Calculate emission species/plant-dependent mapped emission factors and other important coefficients for gamma terms
-        call canopy_biop(EMI_IND, LU_OPT, VTYPE, EF, CT1, CEO, ANEW, AGRO, AMAT, AOLD, ROOTA, ROOTB)
+        call canopy_biop(EMI_IND, LU_OPT, VTYPE, EF, CT1, CEO, ANEW, AGRO, AMAT, AOLD, ROOTA, ROOTB, CAQ, TAQ, DTAQ, &
+            CHT, THT, DTHT, CLT, TLT, DTLT, CHW, THW, DTHW)
 
+!        print*,'CHT=',CHT,'THT=',THT,'DTHT=',DTHT
+!        print*,'CLT=',CLT,'TLT=',TLT,'DTLT=',DTLT
+!        print*,'CHW=',CHW,'THW=',THW,'DTHW=',DTHW
         E_OPT = CEO * EXP(0.05_rk * (TLEAF24_AVE-297.0_rk)) * EXP(0.05_rk * (TLEAF240_AVE-297.0_rk))
 
 ! Calculate gamma (activity) values for average Tleaf (Clifton et al., 2022; based on Guenther et al. 2012)
@@ -205,6 +249,18 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
         GAMMALEAFAGE = GET_GAMMA_LEAFAGE(LEAFAGEOPT, PASTLAI, CURRENTLAI, TSTEPLAI, TABOVECANOPY, ANEW, AGRO, AMAT, AOLD)
         !end do
 
+! Get AQ Stress Factor
+        GAMMAAQ = GET_GAMMA_AQ(AQOPT, W126_REF, W126_SET, CAQ, TAQ, DTAQ)
+
+! Get HT Stress Factor
+        GAMMAHT = GET_GAMMA_HT(HTOPT, DAILY_MAXT2, CHT, THT, DTHT)
+
+! Get LT Stress Factor
+        GAMMALT = GET_GAMMA_LT(LTOPT, DAILY_MINT2, CLT, TLT, DTLT)
+
+! Get HW Stress Factor
+        GAMMAHW = GET_GAMMA_HW(HWOPT, DAILY_MAXWS10, CHW, THW, DTHW)
+
 ! Get canopy loss factor (only used in vertical summing options and empirical formulation and parameters based on isoprene)
 ! Note:  Allowed for other BVOCs but use caution when applying to compare with above canopy flux observations
 
@@ -240,7 +296,7 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
                         FLAI(i) = FLAI(MODLAYS-1)
                     end if
                     EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * GAMMACO2 * CCE * &
-                        GAMMALEAFAGE * GAMMASOIM  ! (ug m-3 hr-1)
+                        GAMMALEAFAGE * GAMMASOIM * GAMMAAQ * GAMMAHT * GAMMALT * GAMMAHW ! (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
@@ -264,7 +320,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
                 end if
             end do
             EMI_OUT(SIZE(ZK)) = LAI * EF * SUM(GammaTLEAF_AVE(1:LAYERS) * GammaPPFD_AVE(1:LAYERS) * &
-                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM * CANLOSS_FAC !put into top model layer (ug m-2 hr-1)
+                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM * &
+                GAMMAAQ * GAMMAHT * GAMMALT * GAMMAHW * CANLOSS_FAC !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
@@ -296,7 +353,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
                 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 * GAMMALEAFAGE * GAMMASOIM * CANLOSS_FAC    !put into top model layer (ug m-2 hr-1)
+                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM * &
+                GAMMAAQ * GAMMAHT * GAMMALT * GAMMAHW * CANLOSS_FAC    !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
@@ -307,7 +365,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, &
                 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 * GAMMALEAFAGE * GAMMASOIM * CANLOSS_FAC    !put into top model layer (ug m-2 hr-1)
+                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM * &
+                GAMMAAQ * GAMMAHT * GAMMALT * GAMMAHW  * CANLOSS_FAC    !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'