SUEWS API Site
Documentation of SUEWS source code
suews_phys_biogenco2.f95
Go to the documentation of this file.
1 module co2_module
2  implicit none
3 
4 contains
5 !========================================================================================
6 ! Created by HCW Aug 2016 to calculate biogenic component of CO2 flux.
7 ! This subroutine is still under development and in the equations there might be bugs and
8 ! the code is not well commented.
9 !
10 ! Last modified:
11 ! MH 18 Feb 2019 - Added new method to calculate photosynthesis with conductance parameters and added option to model with
12 ! 2 meter temperature
13 ! MH 20 Jun 2017 - Tidied and renamed from SUEWS_CO2.f95 to SUEWS_CO2Biogen.f95
14 ! HCW 11 Apr 2017 - Tidied and merged with LJ code
15 ! LJ 6 Apr 2017 - Minimum limit for soil respiration with BiogenCO2Choice = 2 was set to 0.6 umol m-2 s-1
16 ! - Choice for non-rectancular hyperbola to calculate the biogenic CO2 flux added (BiogenCO2Choice = 2)
17 ! (Bellucco et al. 2017, Agric. Forest. Met. 236, 113-122).
18 ! - Both the "local Helsinki model" (BiogenCO2Choice = 2)) and "general model" BiogenCO2Choice = 3) are implemented
19 ! - Snow fraction added to the calculation of active vegetation fraction and the soil respiration
20 !
21 ! To Do:
22 ! - Active vegetation goes to zero with LAI minimum, but this needs to be changed so some minimum value
23 ! especially in the case of evergreentrees
24 ! - Move some of the parameters to input files
25 ! - Now on weekend nighttime population is used throughout the day. Do we need extra column in SiteSelect for daytime weekend population?
26 
27 ! EmissionsMethod:
28 ! 11-16 - Rectangular hyperbola (Ruimy, Schmid, Flanagan)
29 ! 21-26 - Non-rectangular hyperbola, Helsinki (Bellucco et al. 2017)
30 ! 31-36 - Non-rectangular hyperbola, general (Bellucco et al. 2017)
31 ! 41-46 - Rectangular hyperbola, calculated with conductance parameters (Järvi et al., 2019)
32 !========================================================================================
33  SUBROUTINE co2_biogen( &
34  alpha_bioCO2, alpha_enh_bioCO2, avkdn, beta_bioCO2, beta_enh_bioCO2, BSoilSurf, &! input:
35  ConifSurf, DecidSurf, dectime, EmissionsMethod, gfunc, gfunc2, GrassSurf, gsmodel, &
36  id, it, ivConif, ivDecid, ivGrass, LAI_id, LAIMin, LAIMax, min_res_bioCO2, nsurf, &
37  NVegSurf, resp_a, resp_b, sfr, SnowFrac, t2, Temp_C, theta_bioCO2, &
38  Fc_biogen, Fc_photo, Fc_respi)! output:
39 
40  IMPLICIT NONE
41 
42  INTEGER, INTENT(in):: &
43  EmissionsMethod, id, it, ivConif, ivDecid, ivGrass, ConifSurf, DecidSurf, GrassSurf, BSoilSurf, &
44  nsurf, nvegSurf, gsmodel
45 
46  REAL(KIND(1d0)), INTENT(in):: &
47  avkdn, &
48  dectime, &
49  gfunc, &
50  gfunc2, &
51  Temp_C, &
52  t2
53 
54  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in):: &
55  sfr, & !Surface fractions [-]
56  SnowFrac
57 
58  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in):: &
59  LAIMin, LAIMax, & ! [m2 m-2]
60  alpha_bioCO2, &
61  beta_bioCO2, &
62  theta_bioCO2, &
63  alpha_enh_bioCO2, &
64  beta_enh_bioCO2, &
65  resp_a, &
66  resp_b, &
67  min_res_bioCO2, &
68  LAI_id
69 
70  REAL(KIND(1D0)), INTENT(out):: &
71  Fc_biogen, &
72  Fc_respi, Fc_photo
73 
74  INTEGER:: iv ! counter
75 
76  REAL(KIND(1d0)):: &
77  PAR_umolm2s1, &
78  Bellucco2017_Pho, & ! Photosynthesis (Bellucco et al. 2016)
79  Bellucco2017_Res, & ! Respiration (Bellucco et al. 2016)
80  Bellucco2017_Res_surf, &! Respiration for each vegetated surface
81  VegFracSum ! Sum of vegetation fractions without water. Could be moved elsewhere later.
82 
83  REAL(KIND(1d0)), DIMENSION(nvegsurf):: &
84  active_veg_fr, & ! Active vegetation fraction
85  active_veg_fr0, & ! Active vegetation fraction without LAI
86  Fc_photo_surf, & ! Photosynthesis for each vegetated surface
87  Bellucco2017_Pho_surf, & ! Photosynthesis for each vegetated surface
88  alpha_bioCO2_v2, &
89  beta_bioCO2_v2, &
90  theta_bioCO2_v2
91 
92  REAL(KIND(1d0)), PARAMETER :: &
93  JtoumolPAR = 4.6, &
94  kdntopar = 0.46
95 
96  !-----------------------------------------------------------------------
97  ! Calculate PAR from Kdown ----------------------------------
98  par_umolm2s1 = jtoumolpar*kdntopar*avkdn
99  vegfracsum = sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf)
100 
101  ! Calculate active vegetation surface -----------------------
102  !Now this is zero always when LAI in its minimum values. This needs to vary between the summer time maximum and some minimum value
103  !especially in the case of evergreen trees (i.e. in early March LAI can be in its minimum value, but air temperature and radiation
104  !such that uptake can take place)
105  DO iv = ivconif, ivgrass !For vegetated surfaces. Snow included although quite often LAI will be in its minimum when snow on ground
106 
107  active_veg_fr(iv) = (sfr(iv + 2)*(1 - snowfrac(iv + 2)))*(lai_id(iv)/laimax(iv))
108  active_veg_fr0(iv) = (sfr(iv + 2)*(1 - snowfrac(iv + 2)))*lai_id(iv)
109 
110  ENDDO
111 
112  IF (emissionsmethod >= 11 .AND. emissionsmethod <= 16) THEN ! Rectangular hyperbola
113  ! Calculate carbon uptake due to photosynthesis -------------
114  fc_photo = 0
115  DO iv = ivconif, ivgrass
116  fc_photo_surf(iv) = -beta_bioco2(iv)*alpha_bioco2(iv)*par_umolm2s1/(alpha_bioco2(iv)*par_umolm2s1 + beta_bioco2(iv))
117  ! For active vegetation fraction only
118  fc_photo = fc_photo + fc_photo_surf(iv)*active_veg_fr(iv) !umol m-2 s-1
119  ENDDO
120 
121  ELSEIF (emissionsmethod >= 21 .AND. emissionsmethod <= 26) THEN !Local model, Bellucco et al. (2017)
122  ! Calculate carbon uptake due to photosynthesis -------------
123  bellucco2017_pho = 0
124  DO iv = ivconif, ivgrass
125  bellucco2017_pho_surf(iv) = -(1/(2*theta_bioco2(iv)) &
126  *(alpha_bioco2(iv)*par_umolm2s1 + beta_bioco2(iv) &
127  - sqrt((alpha_bioco2(iv)*par_umolm2s1 + beta_bioco2(iv))**2 &
128  - 4*alpha_bioco2(iv)*beta_bioco2(iv)*theta_bioco2(iv)*par_umolm2s1)))
129 
130  ! For active vegetation fraction only
131  bellucco2017_pho = bellucco2017_pho + bellucco2017_pho_surf(iv)*active_veg_fr(iv)
132  ENDDO
133  fc_photo = bellucco2017_pho
134 
135  ELSEIF (emissionsmethod >= 31 .AND. emissionsmethod <= 36) THEN !General model, Bellucco et al. (2017)
136  ! Not currently recommended as includes also some anthropogenic impacts. Should maybe be separate from other biogen choices?
137  ! Alpha and beta calculated as a function of vegetation cover fraction
138 
139  !Different alpha, beta and theta vegetation cover values should be same in BiogenCO2Method = 3.
140  IF (alpha_bioco2(ivconif) == alpha_bioco2(ivdecid) .AND. alpha_bioco2(ivconif) == alpha_bioco2(ivgrass) .AND. &
141  beta_bioco2(ivconif) == beta_bioco2(ivdecid) .AND. beta_bioco2(ivconif) == beta_bioco2(ivgrass) .AND. &
142  theta_bioco2(ivconif) == theta_bioco2(ivdecid) .AND. theta_bioco2(ivconif) == theta_bioco2(ivgrass)) THEN
143 
144  !Because different alpha, beta and theta values are same - only one vegetation type is needed.
145  alpha_bioco2_v2(ivconif) = alpha_bioco2(ivconif) + alpha_enh_bioco2(ivconif)* &
146  (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) + sfr(bsoilsurf)) !umol CO2 umol photons-1
147  beta_bioco2_v2(ivconif) = -beta_bioco2(ivconif) + beta_enh_bioco2(ivconif)* &
148  (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) + sfr(bsoilsurf)) !umol m^-2 s^-1
149 
150  !Photosynthesis
151  bellucco2017_pho = -(1/(2*theta_bioco2(ivconif)) &
152  *(alpha_bioco2_v2(ivconif)*par_umolm2s1 &
153  + beta_bioco2_v2(ivconif) &
154  - sqrt((alpha_bioco2_v2(ivconif)*par_umolm2s1 + beta_bioco2_v2(ivconif))**2 - 4* &
155  alpha_bioco2_v2(ivconif)*beta_bioco2_v2(ivconif)*theta_bioco2(ivconif)*par_umolm2s1)))
156 
157  ELSE !If values are not same, then weighted average is calculated.
158  ! CALL ErrorHint(74, 'Check values in SUEWS_BiogenCO2.txt: ', notUsed, notUsed, notUsedI)
159 
160  ! Weighted averages
161  alpha_bioco2_v2(ivconif) = (alpha_bioco2(ivconif)*sfr(conifsurf)/vegfracsum + &
162  alpha_bioco2(ivdecid)*sfr(decidsurf)/vegfracsum &
163  + alpha_bioco2(ivgrass)*sfr(grasssurf)/vegfracsum) &
164  /(alpha_bioco2(ivconif) + alpha_bioco2(ivdecid) + alpha_bioco2(ivgrass))
165  beta_bioco2_v2(ivconif) = (beta_bioco2(ivconif)*sfr(conifsurf)/vegfracsum + &
166  beta_bioco2(ivdecid)*sfr(decidsurf)/vegfracsum &
167  + beta_bioco2(ivgrass)*sfr(grasssurf)/vegfracsum) &
168  /(beta_bioco2(ivconif) + beta_bioco2(ivdecid) + beta_bioco2(ivgrass))
169  theta_bioco2_v2(ivconif) = (theta_bioco2(ivconif)*sfr(conifsurf)/vegfracsum + &
170  theta_bioco2(ivdecid)*sfr(decidsurf)/vegfracsum &
171  + theta_bioco2(ivgrass)*sfr(grasssurf)/vegfracsum) &
172  /(theta_bioco2(ivconif) + theta_bioco2(ivdecid) + theta_bioco2(ivgrass))
173 
174  ! Using weighted average values to calculate alpha and beta as a function of vegetation cover fraction
175  alpha_bioco2_v2(ivconif) = alpha_bioco2_v2(ivconif) + alpha_enh_bioco2(ivconif)* &
176  (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) + sfr(bsoilsurf)) !umol CO2 umol photons-1
177  beta_bioco2_v2(ivconif) = -beta_bioco2_v2(ivconif) + beta_enh_bioco2(ivconif)* &
178  (sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) + sfr(bsoilsurf)) !umol m^-2 s^-1
179 
180  !Photosynthesis
181  bellucco2017_pho = -(1/(2*theta_bioco2_v2(ivconif))*( &
182  alpha_bioco2_v2(ivconif)*par_umolm2s1 + beta_bioco2_v2(ivconif) &
183  - sqrt((alpha_bioco2_v2(ivconif)*par_umolm2s1 + beta_bioco2_v2(ivconif))**2 &
184  - 4*alpha_bioco2_v2(ivconif)*beta_bioco2_v2(ivconif)*theta_bioco2_v2(ivconif) &
185  *par_umolm2s1)))
186 
187  ENDIF
188  ! Calculate carbon uptake due to photosynthesis -------------
189  fc_photo = bellucco2017_pho*active_veg_fr(conifsurf - 2) + &
190  bellucco2017_pho*active_veg_fr(decidsurf - 2) + &
191  bellucco2017_pho*active_veg_fr(grasssurf - 2)
192 
193  ELSEIF (emissionsmethod >= 41 .AND. emissionsmethod <= 46) THEN ! Conductance parameters
194  !Dependance to incoming shortwave radiation, specific humidity deficit, air temperature and soil moisture deficit
195  fc_photo = 0
196  DO iv = ivconif, ivgrass
197  fc_photo = fc_photo + active_veg_fr0(iv)*beta_bioco2(iv)
198  ENDDO
199 
200  IF (gsmodel == 1 .OR. gsmodel == 2) THEN !With air temperature
201  fc_photo = -fc_photo*gfunc
202  ELSEIF (gsmodel == 3 .OR. gsmodel == 4) THEN !With modelled 2 meter temperature
203  fc_photo = -fc_photo*gfunc2
204  ENDIF
205 
206  ENDIF
207 
208  ! Calculate carbon emissions due to respiration -------------
209  bellucco2017_res = 0.0
210  bellucco2017_res_surf = 0.0
211  IF (vegfracsum > 0.01) THEN
212  DO iv = ivconif, ivgrass
213  IF (sfr(2 + iv) > 0.005) THEN
214  IF (gsmodel == 1 .OR. gsmodel == 2) THEN !With air temperature
215  bellucco2017_res_surf = max(min_res_bioco2(iv), resp_a(iv)*exp(resp_b(iv)*temp_c))
216  ELSEIF (gsmodel == 3 .OR. gsmodel == 4) THEN !With modelled 2 meter temperature
217  bellucco2017_res_surf = max(min_res_bioco2(iv), resp_a(iv)*exp(resp_b(iv)*t2))
218  ENDIF
219  ! Only for vegetation fraction
220  bellucco2017_res = bellucco2017_res + bellucco2017_res_surf*sfr(2 + iv)/vegfracsum
221  ENDIF
222  ENDDO
223  ENDIF
224  fc_respi = bellucco2017_res*(sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) + sfr(bsoilsurf))
225  ! Combine to find biogenic CO2 flux
226  fc_biogen = fc_photo + fc_respi
227 
228  RETURN
229 
230  END SUBROUTINE co2_biogen
231 !========================================================================================
232 
233 end module co2_module
subroutine co2_biogen(alpha_bioCO2, alpha_enh_bioCO2, avkdn, beta_bioCO2, beta_enh_bioCO2, BSoilSurf, ConifSurf, DecidSurf, dectime, EmissionsMethod, gfunc, gfunc2, GrassSurf, gsmodel, id, it, ivConif, ivDecid, ivGrass, LAI_id, LAIMin, LAIMax, min_res_bioCO2, nsurf, NVegSurf, resp_a, resp_b, sfr, SnowFrac, t2, Temp_C, theta_bioCO2, Fc_biogen, Fc_photo, Fc_respi)