Merge pull request noaa-oar-arl#114 from noaa-oar-arl/gamma_sm

Gamma sm

README.md

@@ -149,6 +149,11 @@ The Canopy-App input data in [Table 2](#table-2-canopy-app-required-input-variab
 | `spfh2m`                         | 2-meter specific humidity (kg/kg)           | UFS NOAA/GFSv16                                    |
 | `hpbl`                           | Height of the planetary boundary layer (m)  | UFS NOAA/GFSv16                                    |
 | `prate_ave`                      | Average mass precipitation rate (kg m-2 s-1) | UFS NOAA/GFSv16                                   |
+| `soilw1`                         | Volumetric soil moisture in layer 1 (m3 m-3) | UFS NOAA/GFSv16                                   |
+| `soilw2`                         | Volumetric soil moisture in layer 2 (m3 m-3) | UFS NOAA/GFSv16                                   |
+| `soilw3`                         | Volumetric soil moisture in layer 3 (m3 m-3) | UFS NOAA/GFSv16                                   |
+| `soilw4`                         | Volumetric soil moisture in layer 4 (m3 m-3) | UFS NOAA/GFSv16                                   |
+| `wilt`                           | Wilting point (proportion)                  | UFS NOAA/GFSv16                                    |
 | **External Canopy Variables**    | **Variable Description and Units**          | **Data Source/Reference (if necessary)**           |
 | `ch`                             | Canopy height (m)                    | Globally extended GEDI data. Data Period=2020. Data frequency=Annual. ([Lang et al., 2023](https://doi.org/10.1038/s41559-023-02206-6)) |
 | `clu`                            | Canopy clumping index (dimensionless)       | GriddingMachine/MODIS. Data Period=2001-2017 Climatology. Data frequency=Monthly. ([Wei et al., 2019](https://doi.org/10.1016/j.rse.2019.111296)). Extended globally for high latitudes using methods described [here](https://gmuedu-my.sharepoint.com/:w:/g/personal/whung_gmu_edu/EdglXmW2kzBDtDj1xV0alGcB1Yo2I8hzdyWGVGB2YOTfgw). |
@@ -249,6 +254,8 @@ You can also [generate global inputs using Python (see python/global_data_proces
 | `co2_set`       | user-set real value of atmospheric co2 concentration (ppmv) (only used if `co2_opt=0` or `co2_opt=1`) |
 | `leafage_opt`   | user-set options for applying leaf-age response to biogenic VOC emissions based on [Guenther et al. 2006](https://doi.org/10.5194/acp-6-3181-2006) (default is off i.e., `leafage_opt=1`, the corresponding $\gamma$ is set to 1). If turned on (`leafage_opt=0`), leafage $\gamma$ is calculated and the lai_tstep option needs to be set to ensure correct interpolation in this leafage_opt calculation. |
 | `lai_tstep`     | user-defined options for the number of seconds in the interval at which LAI (Leaf Area Index) data is provided to the model. For instance, if LAI data is given on a daily basis, lai_tstep would be set to the number of seconds in a day (86,400 seconds). If LAI data is provided monthly, then lai_tstep would represent the total number of seconds in that month (e.g., 2,592,000 seconds for a 30-day month). This parameter helps in determining the frequency of LAI input and is crucial for interpolating LAI values to the model's hourly timesteps when the model's timestep (time_intvl) is smaller than the LAI input interval.  |
+| `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`. |
 
 **\*\*** 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

@@ -73,4 +73,9 @@
   co2_set     =  400.0
   lai_tstep   =  86400
   leafage_opt =  1
+  soim_opt    =  1
+  soild1      =  5.0
+  soild2      =  25.0
+  soild3      =  70.0
+  soild4      =  150.0
 /

input/point_file_20220630.sfcf023.txt

@@ -1,2 +1,2 @@
-lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave
-34.03,272.11,20.8692,-0.7485,0.3286,0.5214,3.7059,4,0.5124,0.1107,0.0416,0.8260,12.2366,0.0000,10.00,3,99731.6094,10.1607,-9.8159,292.9766,293.2202,0.0146,14.0288,0.0001
+lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave,soilw1,soilw2,soilw3,soilw4,wilt
+34.03,272.11,20.8692,-0.7485,0.3286,0.5214,3.7059,4,0.5124,0.1107,0.0416,0.8260,12.2366,0.0000,10.00,3,99731.6094,10.1607,-9.8159,292.9766,293.2202,0.0146,14.0288,0.0001,0.2591,0.2053,0.2134,0.2132,0.0470

input/point_file_20220701.sfcf000.txt

@@ -1,2 +1,2 @@
-lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave
-34.03,272.11,20.8692,0.1580,0.8375,0.5105,3.6504,4,0.5124,0.1548,0.2357,0.8260,-181.1201,0.0000,10.00,3,99809.6719,122.7623,1.8252,294.9718,294.9180,0.0157,85.5025,0.0000
+lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave,soilw1,soilw2,soilw3,soilw4,wilt
+34.03,272.11,20.8692,0.1580,0.8375,0.5105,3.6504,4,0.5124,0.1548,0.2357,0.8260,-181.1201,0.0000,10.00,3,99809.6719,122.7623,1.8252,294.9718,294.9180,0.0157,85.5025,0.0000,0.2487,0.2256,0.2283,0.2307,0.0470

input/point_file_20220701.sfcf001.txt

@@ -1,2 +1,2 @@
-lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave
-34.03,272.11,20.8692,-0.0656,0.8462,0.5105,3.6504,4,0.5124,0.1866,0.4319,0.8260,-16.3281,0.0000,10.00,3,99854.1172,336.8927,35.6389,296.9243,296.5100,0.0162,234.9751,0.0000
+lat,lon,ch,ugrd10m,vgrd10m,clu,lai,vtype,canfrac,fricv,csz,sfcr,mol,frp,href,sotyp,pressfc,dswrf,shtfl,tmpsfc,tmp2m,spfh2m,hpbl,prate_ave,soilw1,soilw2,soilw3,soilw4,wilt
+34.03,272.11,20.8692,-0.0656,0.8462,0.5105,3.6504,4,0.5124,0.1866,0.4319,0.8260,-16.3281,0.0000,10.00,3,99854.1172,336.8927,35.6389,296.9243,296.5100,0.0162,234.9751,0.0000,0.2477,0.2259,0.2283,0.2307,0.0470

src/canopy_bioemi_mod.F90

@@ -9,6 +9,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
         PPFD_SHADE, TLEAF_SUN, TLEAF_SHADE, TLEAF_AVE, TEMP2, LU_OPT, &
         VTYPE, MODRES, CCE, VERT, CO2OPT, CO2SET, &
         LEAFAGEOPT, PASTLAI, CURRENTLAI, TSTEPLAI, &
+        SOIMOPT, SOIM1, SOIM2, SOIM3, SOIM4, SOID1, SOID2, SOID3, &
+        SOID4, WILT, &
         MODLAYS, EMI_IND, EMI_OUT)
 
 !-----------------------------------------------------------------------
@@ -38,8 +40,10 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
 !     BVOCs
 !-----------------------------------------------------------------------
 !-----------------------------------------------------------------------
-        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
+        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
         use canopy_bioparm_mod
         use canopy_tleaf_mod
 
@@ -64,6 +68,16 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
         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 )       :: SOIMOPT         ! Option for soil moisture factor
+        REAL(RK),    INTENT( IN )       :: SOIM1           ! Volumetric soil moisture layer 1 [m3/m3]
+        REAL(RK),    INTENT( IN )       :: SOIM2           ! Volumetric soil moisture layer 2 [m3/m3]
+        REAL(RK),    INTENT( IN )       :: SOIM3           ! Volumetric soil moisture layer 3 [m3/m3]
+        REAL(RK),    INTENT( IN )       :: SOIM4           ! Volumetric soil moisture layer 4 [m3/m3]
+        REAL(RK),    INTENT( IN )       :: SOID1           ! Soil depth layer 1 [cm]
+        REAL(RK),    INTENT( IN )       :: SOID2           ! Soil depth layer 2 [cm]
+        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 )       :: LEAFAGEOPT      ! leafage_opt (0= ON, 1= off i.e. GAMMALEAFAGE =1, in canopy_readnml.F90)
         REAL(RK),    INTENT( IN )       :: PASTLAI           ! Past LAI [cm2/cm2]
@@ -113,6 +127,11 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
         REAL(RK) :: AMAT
         REAL(RK) :: AOLD
 
+        ! 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) :: GAMMACO2                       ! CO2 inhibition factor (isoprene only)
 
         REAL(RK) :: GAMMALEAFAGE !(SIZE(ZK))                 ! LEAF AGE factor
@@ -129,8 +148,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
         TLEAF24_AVE    = TLEAF_AVE
         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
-        call canopy_biop(EMI_IND, LU_OPT, VTYPE, EF, CT1, CEO, ANEW, AGRO, AMAT, AOLD)  !Update: Added leaf age empirical coeff's
+! 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)
 
         E_OPT = CEO * EXP(0.05_rk * (TLEAF24_AVE-297.0_rk)) * EXP(0.05_rk * (TLEAF240_AVE-297.0_rk))
 
@@ -168,6 +187,10 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
             GAMMACO2 = 1.0_rk
         end if
 
+! Get Soil Moisture Factor
+        GAMMASOIM =  GET_GAMMA_SOIM(SOIMOPT,SOIM1,SOIM2,SOIM3,SOIM4,SOID1,SOID2,SOID3,SOID4,WILT, &
+            ROOTA,ROOTB)
+
 ! Get LEAF AGE factor
         TABOVECANOPY  = TEMP2   !TEMP2 (above air temp) for TABOVECANOPY
         !do i=1, SIZE(ZK)
@@ -186,7 +209,8 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
                     else
                         FLAI(i) = FLAI(MODLAYS-1)
                     end if
-                    EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * GAMMACO2 * CCE * GAMMALEAFAGE  ! (ug m-3 hr-1)
+                    EMI_OUT(i) = FLAI(i) * EF * GammaTLEAF_AVE(i) * GammaPPFD_AVE(i) * GAMMACO2 * CCE * &
+                        GAMMALEAFAGE * GAMMASOIM  ! (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
@@ -210,7 +234,7 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
                 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 !put into top model layer (ug m-2 hr-1)
+                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM !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
@@ -242,7 +266,7 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
                 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    !put into top model layer (ug m-2 hr-1)
+                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM    !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
@@ -253,7 +277,7 @@ SUBROUTINE CANOPY_BIO( ZK, FCLAI, FCH, LAI, FSUN, PPFD_SUN, &
                 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    !put into top model layer (ug m-2 hr-1)
+                VPGWT(1:LAYERS)) * GAMMACO2 * CCE * GAMMALEAFAGE * GAMMASOIM    !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'