SUEWS API Site
Documentation of SUEWS source code
suews_phys_anemsn.f95
Go to the documentation of this file.
2  implicit none
3 
4  !===================================================================================
5  !Simple Anthropogenic Heat and Carbon Dioxide Parameterization routines
6  !This subroutine is still under development and in the equations concerning CO2 fluxes
7  !there are bugs and missing comments.
8  !Last modified
9  ! TS 01 Jul 2018 - replaced annual HDD array with a simplied daily array HDD_id_use
10  ! MH 29 Jun 2017 - Finalised the code to calculate the anthropogenic emissions of heat and CO2
11  ! HCW 21 Apr 2017 - renamed from SUEWS_SAHP.f95. Now includes CO2 fluxes as well as QF.
12  ! Tidied code. Combined three subroutines into 1 to avoid repetition.
13  ! HCW 24 Feb 2017 - Added new anthropogenic heat flux calculation (AnthropEmissionsMethod=3)
14  ! HCW 25 Aug 2016 - Outputs base QF (part without temp. dependence)
15  ! LJ 27 Jan 2016 - Removal of Tabs
16  ! HCW 20 Jan 2015 - v2015 applies a profile at each model timestep
17  ! these have been interpolated from the hourly profile input data (SUEWS_Profiles)
18  ! vthey are now normalised (sum to 1) in InitializeSurfaceCharacteristics
19  ! N.B. previous versions were not applying weekday/weekend profiles correctly
20  !
21  ! EmissionsMethod = 0 - Use values in met forcing file, or default QF
22  ! EmissionsMethod = 1 - Method according to Loridan et al. (2011) : SAHP
23  ! EmissionsMethod = 2 - Method according to Jarvi et al. (2011) : SAHP_2
24  !
25  !===================================================================================
26 
27 contains
28 
29  SUBROUTINE anthropogenicemissions( &
30  CO2PointSource, EmissionsMethod, &
31  it, imin, DLS, DayofWeek_id, &
32  EF_umolCO2perJ, FcEF_v_kgkm, EnEF_v_Jkm, TrafficUnits, &
33  FrFossilFuel_Heat, FrFossilFuel_NonHeat, &
34  MinFCMetab, MaxFCMetab, MinQFMetab, MaxQFMetab, &
35  PopDensDaytime, PopDensNighttime, &
36  Temp_C, HDD_id, Qf_A, Qf_B, Qf_C, &
37  AH_MIN, AH_SLOPE_Heating, AH_SLOPE_Cooling, &
38  BaseT_Heating, BaseT_Cooling, &
39  TrafficRate, &
40  QF0_BEU, QF_SAHP, &
41  Fc_anthro, Fc_metab, Fc_traff, Fc_build, Fc_point, &
42  AHProf_24hr, HumActivity_24hr, TraffProf_24hr, PopProf_24hr, SurfaceArea)
43 
44  ! Simple anthropogenic heat parameterisation and co2 calculation
45  ! Calculates QF_SAHP and Fc_anthro
46 
47  IMPLICIT NONE
48 
49  INTEGER, INTENT(in):: &
50  EmissionsMethod, &
51  it, & !Hour
52  imin, & !Minutes
53  DLS !day lightsavings =1+1h=0
54 
55  INTEGER, DIMENSION(3), INTENT(in)::DayofWeek_id !1 - day of week; 2 - month; 3 - season
56 
57  REAL(KIND(1d0)), DIMENSION(12), INTENT(in):: HDD_id !Heating Degree Days (see SUEWS_DailyState.f95)
58 
59  REAL(KIND(1d0)), DIMENSION(2), INTENT(in):: &
60  Qf_A, Qf_B, Qf_C, & !Qf coefficients
61  AH_MIN, & !Minimum anthropogenic heat flux
62  AH_SLOPE_Heating, & !Slope of the anthropogenic heat flux calculation
63  AH_SLOPE_Cooling, &
64  FcEF_v_kgkm, & !CO2 Emission factor
65  PopDensDaytime, & !Daytime population density [ha-1] (i.e. workers)
66  BaseT_Heating, & !Critical temperature
67  BaseT_Cooling, & !Critical cooling temperature
68  TrafficRate, & !Traffic rate
69  QF0_BEU
70 
71  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in):: AHProf_24hr
72  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in):: HumActivity_24hr
73  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in):: TraffProf_24hr
74  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in):: PopProf_24hr
75 
76  REAL(KIND(1D0)), INTENT(in):: &
77  CO2PointSource, & !Point source [kgC day-1]
78  EF_umolCO2perJ, &
79  EnEF_v_Jkm, & !Energy Emission factor
80  TrafficUnits, & !Traffic Units choice
81  FrFossilFuel_Heat, & !Fraction of Fossil Fuel heat
82  FrFossilFuel_NonHeat, &!Fraction of Fossil Fuel non heat
83  MinFCMetab, & !Minimum FC Metabolism
84  MaxFCMetab, & !Maximum FC Metabolism
85  MinQFMetab, & !Minimum QF Metabolism
86  MaxQFMetab, & !Maximum QF Metabolism
87  PopDensNighttime, & !Nighttime population density [ha-1] (i.e. residents)
88  SurfaceArea, & !Surface area [m-2]
89  Temp_C !Air temperature
90 
91  REAL(KIND(1D0)), INTENT(out):: &
92  QF_SAHP, &
93  Fc_anthro, &
94  Fc_metab, Fc_traff, Fc_build, Fc_point
95 
96  INTEGER:: &
97  iu, & !1=weekday OR 2=weekend
98  ih
99 
100  REAL(KIND(1D0)):: &
101  get_Prof_SpecTime_inst, &
102  get_Prof_SpecTime_mean, & !external function to get profile value at specified timestamp
103  HDD_daily, & !daily HDD
104  CDD_daily, & !daily HDD
105  Tair_avg_daily, & !daily mean air temperature
106  NumCapita(2), &
107  DP_x_RhoPop, DP_x_RhoPop_traff, &
108  QF_build, QF_metab, QF_traff, &
109  QF_SAHP_base, & !Anthropogenic heat flux calculated by SAHP (temp independent part)
110  QF_SAHP_heating, & !Anthropogenic heat flux calculated by SAHP (heating part only)
111  QF_SAHP_cooling, & !AC contribution
112  PopDorNorT, & ! Population
113  ActDorNorT, & ! Human activity
114  TraffDorNorT, & ! Traffic
115  AHDorNorT ! Anthropogenic heat
116 
117  ! initialisation
118  qf_traff = 0
119  qf_sahp_heating = 0
120  qf_sahp_cooling = 0
121 
122  ! Transfer HDD/CDD values to local explict variables
123  hdd_daily = hdd_id(7)
124  cdd_daily = hdd_id(8)
125  ! NB: temporarily use 5-day running mean to test the performance
126  tair_avg_daily = hdd_id(10)
127 
128  ! Tair_avg_daily= HDD_id_use(3) ! this is daily
129 
130  ! use average of both daytime and nighttime population values
131  ! TS 20 Sep 2019: moved from here to simplify the interface
132  IF (popdensdaytime(1) >= 0 .AND. popdensnighttime >= 0) numcapita(1) = (popdensdaytime(1) + popdensnighttime)/2 !If both, use average
133  IF (popdensdaytime(2) >= 0 .AND. popdensnighttime >= 0) numcapita(2) = (popdensdaytime(2) + popdensnighttime)/2 !If both, use average
134 
135  !NumCapita = (PopDensDaytime + PopDensNighttime)/2
136 
137  !-----------------------------------------------------------------------
138  ! Account for Daylight saving
139  ih = it - dls
140  IF (ih < 0) ih = 23
141 
142  ! Set weekday/weekend counter
143  iu = 1 !Set to 1=weekday
144  IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) iu = 2 !Set to 2=weekend
145 
146  ! Calculate energy emissions and CO2 from human metabolism -------------
147  ! Pop dens (cap ha-1 -> cap m-2) x activity level (W cap-1)
148  ! PopDorNorT = PopProf_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu) ! 1=night, 2=day, 1-2=transition
149  ! ActDorNorT = HumActivity_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu) ! 1=night, 2=day, 1-2=transition
150  ! TraffDorNorT = TraffProf_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu) ! normalise so the AVERAGE of the multipliers is equal to 1
151  ! AHDorNorT = AHProf_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu) ! normalise so the AVERAGE of the multipliers is equal to 1
152  ! PRINT*, ''
153  ! PRINT*, 'AHDorNorT old:',AHDorNorT
154  popdornort = get_prof_spectime_inst(ih, imin, 0, popprof_24hr(:, iu)) ! 1=night, 2=day, 1-2=transition
155  actdornort = get_prof_spectime_inst(ih, imin, 0, humactivity_24hr(:, iu)) ! 1=night, 2=day, 1-2=transition
156  traffdornort = get_prof_spectime_mean(ih, imin, 0, traffprof_24hr(:, iu)) ! normalise so the AVERAGE of the multipliers is equal to 1
157  ahdornort = get_prof_spectime_mean(ih, imin, 0, ahprof_24hr(:, iu)) ! normalise so the AVERAGE of the multipliers is equal to 1
158 
159  ! Diurnal profile times population density [cap ha-1]
160  dp_x_rhopop = ahdornort*numcapita(iu)
161 
162  qf_metab = (popdensnighttime*minqfmetab*((2 - actdornort) + (2 - popdornort))/2 + &
163  popdensdaytime(iu)*maxqfmetab*((actdornort - 1) + (popdornort - 1))/2)/10000 !W m-2
164 
165  fc_metab = (popdensnighttime*minfcmetab*((2 - actdornort) + (2 - popdornort))/2 + &
166  popdensdaytime(iu)*maxfcmetab*((actdornort - 1) + (popdornort - 1))/2)/10000 !umol m-2 s-1
167 
168  !----------------------------------------------------------------------------------------------------
169  !Different methods to calculate the anthopogenic heat and CO2 emissions.
170  !1-3: CO2 emission is calculated from QF which can be calculated with three methods
171  !41-43: CO2 emission is calculated with local information. QF methods are used for housing and human metabolism
172 
173  IF (emissionsmethod == 1 .OR. emissionsmethod == 4 .OR. &
174  emissionsmethod == 11 .OR. emissionsmethod == 14 .OR. &
175  emissionsmethod == 21 .OR. emissionsmethod == 24 .OR. &
176  emissionsmethod == 31 .OR. emissionsmethod == 34 .OR. &
177  emissionsmethod == 41 .OR. emissionsmethod == 44) THEN ! (formerly SAHP_1 subroutine)
178  ! Loridan et al. (2011) JAMC Eq 13: linear relation with air temperature
179  ! Weekday/weekend differences due to profile only
180  ! Now scales with population density
181 
182  IF (temp_c < baset_heating(iu)) THEN
183  ! QF_SAHP = (AH_MIN(iu) + AH_SLOPE_Heating(iu)*(BaseT_Heating(iu) - Temp_C))*AHDorNorT
184  qf_sahp_heating = (ah_slope_heating(iu)*(baset_heating(iu) - temp_c))*ahdornort
185  ELSE
186  qf_sahp_heating = 0
187  ! QF_SAHP = AH_MIN(iu)*AHDorNorT
188  ENDIF
189 
190  ! Need to be checked later, not recommended to use
191  qf_sahp_base = ah_min(iu)*ahdornort ! Temperature-independent contribution
192  qf_sahp_cooling = 0 ! No AC contribution with this method
193 
194  ELSEIF (emissionsmethod == 2 .OR. emissionsmethod == 5 .OR. &
195  emissionsmethod == 12 .OR. emissionsmethod == 15 .OR. &
196  emissionsmethod == 22 .OR. emissionsmethod == 25 .OR. &
197  emissionsmethod == 32 .OR. emissionsmethod == 35 .OR. &
198  emissionsmethod == 42 .OR. emissionsmethod == 45) THEN ! (formerly SAHP_2 subroutine)
199  ! Jarvi et al. (2011) JH Eq 3 using HDD and CDD
200  ! Weekday/weekend differences due to profile and coefficients QF_a,b,c
201  ! Scales with population density
202  ! QF_SAHP = (Qf_a(iu) + Qf_b(iu)*CDD_daily + Qf_c(iu)*HDD_daily)*DP_x_RhoPop !This contains QF from all three sources: buildings, metabolism and traffic!
203  qf_sahp_base = (qf_a(iu))*dp_x_rhopop ! Temperature-independent contribution from buildings, traffic and human metabolism
204  qf_sahp_heating = (qf_c(iu)*hdd_daily)*dp_x_rhopop ! Heating contribution
205  qf_sahp_cooling = (qf_b(iu)*cdd_daily)*dp_x_rhopop ! Cooling (AC) contribution
206 
207  ELSEIF (emissionsmethod == 3 .OR. emissionsmethod == 6 .OR. &
208  emissionsmethod == 13 .OR. emissionsmethod == 16 .OR. &
209  emissionsmethod == 23 .OR. emissionsmethod == 26 .OR. &
210  emissionsmethod == 33 .OR. emissionsmethod == 36 .OR. &
211  emissionsmethod == 43 .OR. emissionsmethod == 46) THEN
212  ! Updated Loridan et al. (2011) method using daily (not instantaneous) air temperature (HDD_id_use(3))
213  ! Linear relation with air temperature
214  ! Weekday/weekend differences due to profile only
215  ! Scales with population density
216 
217  ! Need to be checked later, not recommended to use
218  qf_sahp_base = ah_min(iu)*ahdornort ! Temperature-independent contribution
219 
220  IF (tair_avg_daily < baset_heating(iu)) THEN ! Heating
221  qf_sahp_heating = (ah_slope_heating(iu)*(baset_heating(iu) - tair_avg_daily))*ahdornort ! Heating contribution
222  qf_sahp_cooling = 0
223 
224  ELSEIF (tair_avg_daily > baset_cooling(iu)) THEN ! Air-conditioning
225  qf_sahp_heating = 0
226  qf_sahp_cooling = (ah_slope_cooling(iu)*(tair_avg_daily - baset_cooling(iu)))*ahdornort ! AC contribution
227 
228  ELSE
229  qf_sahp_heating = 0
230  qf_sahp_cooling = 0
231  ENDIF
232 
233  ENDIF
234 
235  ! total QF as the sum of components
236  qf_sahp = qf_sahp_base + qf_sahp_heating + qf_sahp_cooling
237 
238  IF (emissionsmethod >= 1 .AND. emissionsmethod <= 3 .OR. &
239  emissionsmethod >= 11 .AND. emissionsmethod <= 13 .OR. &
240  emissionsmethod >= 21 .AND. emissionsmethod <= 23 .OR. &
241  emissionsmethod >= 31 .AND. emissionsmethod <= 33 .OR. &
242  emissionsmethod >= 41 .AND. emissionsmethod <= 43) THEN
243 
244  ! Calculate QF from buildings.
245  ! First remove (if possibe) human metabolism from the total value given by SAHP.
246  IF ((qf_sahp_base - qf_metab) > 0) THEN
247  qf_build = qf_sahp_base*qf0_beu(iu) + qf_sahp_heating + qf_sahp_cooling !QF0_BEU = QF0_BuildingEnergyUse = Fraction of base value coming from buildings
248  !relative to traffic as metabolism is separately calculated
249  ELSE
250  CALL errorhint(69, 'QF metab exceeds base QF.', qf_metab, qf_sahp_base)
251 
252  qf_build = qf_sahp_heating + qf_sahp_cooling + (qf_sahp_base - qf_metab) !If human metabolism greater than Base QF, remove this from the heating/cooling contribution also
253  ENDIF
254 
255  ! Consider the various components of QF_build to calculate Fc_build
256  ! Assume all A/C electric, so QF_SAHP_ac is not associated with any local CO2 emissions
257  ! HDD part is building energy use, split between electricity (no local emissions CO2) and combustion (CO2) heating...
258  fc_build = qf_sahp_heating*frfossilfuel_heat*ef_umolco2perj
259 
260  ! ... and there is also a temperature-independent contribution from building energy use.
261  IF ((qf_sahp_base - qf_metab) > 0) THEN
262  fc_build = fc_build + qf_sahp_base*qf0_beu(iu)*frfossilfuel_nonheat*ef_umolco2perj
263  ENDIF
264 
265  ! Remainder of temperature-independent part must be traffic
266  !QF_traff = (QF_SAHP_base-QF_metab)*(1.0-QF0_BEU(iu))
267  qf_traff = qf_sahp_base*(1.0 - qf0_beu(iu)) - qf_metab
268 
269  ! Divide QF by energy emission factor and multiply by CO2 factor
270  fc_traff = qf_traff/enef_v_jkm*fcef_v_kgkm(iu)*1e3*1e6/44
271 
272  !CO2 point source [kg C day-1] * 1e3 g kg-1 * 1e6 umol mol-1 /(12 g mol-1 * m-2)
273  IF (co2pointsource > 0) THEN
274  fc_point = co2pointsource*1e3*1e6/(12*60*60*24*surfacearea)
275  ELSE
276  fc_point = 0
277  ENDIF
278 
279  ! Sum components to give anthropogenic CO2 flux [umol m-2 s-1]
280  fc_anthro = fc_metab + fc_traff + fc_build + fc_point
281 
282  ELSEIF (emissionsmethod >= 4 .AND. emissionsmethod <= 6 .OR. &
283  emissionsmethod >= 14 .AND. emissionsmethod <= 16 .OR. &
284  emissionsmethod >= 24 .AND. emissionsmethod <= 26 .OR. &
285  emissionsmethod >= 34 .AND. emissionsmethod <= 36 .OR. &
286  emissionsmethod >= 44 .AND. emissionsmethod <= 46) THEN
287  ! Calculate QF and Fc using building energy use and transport statistics
288 
289  ! Calculate energy released from traffic
290  IF (trafficunits == 1) THEN ! [veh km m-2 day-1]
291  ! Calculate using mean traffic rate [veh km m-2 day-1] * emission factor [J km-1]
292  qf_traff = trafficrate(iu)/(60*60*24)*enef_v_jkm*traffdornort
293  ! Use mean traffic rate [veh km m-2 s-1] * emission factor [kg km-1] * 1e3 g kg-1 /44 g mol-1 * 1e6 umol mol-1
294  fc_traff = trafficrate(iu)/(60*60*24)*fcef_v_kgkm(iu)*1e3*1e6/44*traffdornort
295 
296  ELSEIF (trafficunits == 2) THEN ! [veh km cap-1 day-1]
297  dp_x_rhopop_traff = traffdornort*numcapita(iu)/10000 ! [cap m-2]
298  ! Use mean traffic rate [veh km cap-1 day-1] * emission factor [J km-1]
299  qf_traff = trafficrate(iu)/(60*60*24)*enef_v_jkm*dp_x_rhopop_traff
300  ! Use mean traffic rate [veh km cap-1 day-1] * emission factor [kg km-1] * 1e3 g kg-1 /44 g mol-1 * 1e6 umol mol-1
301  fc_traff = trafficrate(iu)/(60*60*24)*fcef_v_kgkm(iu)*1e3*1e6/44*dp_x_rhopop_traff
302 
303  ELSE ! If TrafficUnits doesn't match possible units
304  CALL errorhint(75, 'Check TrafficUnits', trafficunits)
305 
306  ENDIF
307 
308  ! Energy released from buildings only
309  ! Should buildings have their own profile? Now using population profile
310  qf_build = ((qf_sahp_base*qf0_beu(iu) + qf_sahp_heating + qf_sahp_cooling)/dp_x_rhopop)* &
311  (popdensnighttime*(2 - popdornort) + popdensdaytime(iu)*(popdornort - 1))
312 
313  ! Consider the various components of QF_build to calculate Fc_build
314  fc_build = qf_sahp_heating*frfossilfuel_heat*ef_umolco2perj
315  ! ... and there is also a temperature-independent contribution from building energy use
316  fc_build = fc_build + qf_sahp_base*qf0_beu(iu)*frfossilfuel_nonheat*ef_umolco2perj
317 
318  !CO2 point source [kg C day-1] * 1e3 g kg-1 * 1e6 umol mol-1 /(12 g mol-1 * m-2)
319  IF (co2pointsource > 0) THEN
320  fc_point = co2pointsource*1e3*1e6/(12*60*60*24*surfacearea)
321  ELSE
322  fc_point = 0
323  ENDIF
324 
325  ! Add other QF components to QF_SAHP_base (because if AnthropEmissionsMethod = 4, QF calculated above is building part only)
326  qf_sahp_base = qf_sahp_base + qf_traff + qf_metab
327 
328  ! Sum components to give anthropogenic heat flux [W m-2]
329  qf_sahp = qf_metab + qf_traff + qf_build
330 
331  ! Sum components to give anthropogenic CO2 flux [umol m-2 s-1]
332  fc_anthro = fc_metab + fc_traff + fc_build + fc_point
333 
334  ENDIF
335 
336  RETURN
337  END SUBROUTINE anthropogenicemissions
338 
339 end module anemsn_module
340 !========================================================================================
subroutine anthropogenicemissions(CO2PointSource, EmissionsMethod, it, imin, DLS, DayofWeek_id, EF_umolCO2perJ, FcEF_v_kgkm, EnEF_v_Jkm, TrafficUnits, FrFossilFuel_Heat, FrFossilFuel_NonHeat, MinFCMetab, MaxFCMetab, MinQFMetab, MaxQFMetab, PopDensDaytime, PopDensNighttime, Temp_C, HDD_id, Qf_A, Qf_B, Qf_C, AH_MIN, AH_SLOPE_Heating, AH_SLOPE_Cooling, BaseT_Heating, BaseT_Cooling, TrafficRate, QF0_BEU, QF_SAHP, Fc_anthro, Fc_metab, Fc_traff, Fc_build, Fc_point, AHProf_24hr, HumActivity_24hr, TraffProf_24hr, PopProf_24hr, SurfaceArea)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)