suews_phys_anemsn.f95

Go to the documentation of this file.

MODULE anemsn_module
   IMPLICIT NONE
 
   !===================================================================================
   !Simple Anthropogenic Heat and Carbon Dioxide Parameterization routines
   !This subroutine is still under development and in the equations concerning CO2 fluxes
   !there are bugs and missing comments.
   !Last modified
   ! TS 01 Jul 2018 - replaced annual HDD array with a simplied daily array HDD_id_use
   ! MH 29 Jun 2017 -  Finalised the code to calculate the anthropogenic emissions of heat and CO2
   ! HCW 21 Apr 2017 - renamed from SUEWS_SAHP.f95. Now includes CO2 fluxes as well as QF.
   !                   Tidied code. Combined three subroutines into 1 to avoid repetition.
   ! HCW 24 Feb 2017 - Added new anthropogenic heat flux calculation (AnthropEmissionsMethod=3)
   ! HCW 25 Aug 2016 - Outputs base QF (part without temp. dependence)
   ! LJ 27 Jan 2016  - Removal of Tabs
   ! HCW 20 Jan 2015 - v2015 applies a profile at each model timestep
   !                   these have been interpolated from the hourly profile input data (SUEWS_Profiles)
   !                   vthey are now normalised (sum to 1) in InitializeSurfaceCharacteristics
   !                   N.B. previous versions were not applying weekday/weekend profiles correctly
   !
   ! EmissionsMethod = 0 - Use values in met forcing file, or default QF
   ! EmissionsMethod = 1 - Method according to Loridan et al. (2011) : SAHP
   ! EmissionsMethod = 2 - Method according to Jarvi et al. (2011)   : SAHP_2
   !
   !===================================================================================
 
CONTAINS
 
   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)
 
      ! Simple anthropogenic heat parameterisation and co2 calculation
      ! Calculates QF_SAHP and Fc_anthro
 
      IMPLICIT NONE
 
      INTEGER, INTENT(in) :: &
         EmissionsMethod, &
         it, & !Hour
         imin, & !Minutes
         DLS !day lightsavings =1+1h=0
 
      INTEGER, DIMENSION(3), INTENT(in) :: DayofWeek_id !1 - day of week; 2 - month; 3 - season
 
      REAL(KIND(1D0)), DIMENSION(12), INTENT(in) :: HDD_id !Heating Degree Days (see SUEWS_DailyState.f95)
 
      REAL(KIND(1D0)), DIMENSION(2), INTENT(in) :: &
         Qf_A, Qf_B, Qf_C, & !Qf coefficients
         AH_MIN, & !Minimum anthropogenic heat flux
         AH_SLOPE_Heating, & !Slope of the anthropogenic heat flux calculation
         AH_SLOPE_Cooling, &
         FcEF_v_kgkm, & !CO2 Emission factor
         PopDensDaytime, & !Daytime population density [ha-1] (i.e. workers)
         BaseT_Heating, & !Critical temperature
         BaseT_Cooling, & !Critical cooling temperature
         TrafficRate, & !Traffic rate
         QF0_BEU
 
      REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: AHProf_24hr
      REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: HumActivity_24hr
      REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: TraffProf_24hr
      REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: PopProf_24hr
 
      REAL(KIND(1D0)), INTENT(in) :: &
         CO2PointSource, & !Point source [kgC day-1]
         EF_umolCO2perJ, &
         EnEF_v_Jkm, & !Energy Emission factor
         TrafficUnits, & !Traffic Units choice
         FrFossilFuel_Heat, & !Fraction of Fossil Fuel heat
         FrFossilFuel_NonHeat, & !Fraction of Fossil Fuel non heat
         MinFCMetab, & !Minimum FC Metabolism
         MaxFCMetab, & !Maximum FC Metabolism
         MinQFMetab, & !Minimum QF Metabolism
         MaxQFMetab, & !Maximum QF Metabolism
         PopDensNighttime, & !Nighttime population density [ha-1] (i.e. residents)
         SurfaceArea, & !Surface area [m-2]
         Temp_C !Air temperature
 
      REAL(KIND(1D0)), INTENT(out) :: &
         QF_SAHP, &
         Fc_anthro, &
         Fc_metab, Fc_traff, Fc_build, Fc_point
 
      INTEGER :: &
         iu, & !1=weekday OR 2=weekend
         ih
 
      REAL(KIND(1D0)) :: &
         get_Prof_SpecTime_inst, &
         get_Prof_SpecTime_mean, & !external function to get profile value at specified timestamp
         HDD_daily, & !daily HDD
         CDD_daily, & !daily HDD
         Tair_avg_daily, & !daily mean air temperature
         NumCapita(2), &
         DP_x_RhoPop, DP_x_RhoPop_traff, &
         QF_build, QF_metab, QF_traff, &
         QF_SAHP_base, & !Anthropogenic heat flux calculated by SAHP (temp independent part)
         QF_SAHP_heating, & !Anthropogenic heat flux calculated by SAHP (heating part only)
         QF_SAHP_cooling, & !AC contribution
         PopDorNorT, & ! Population
         ActDorNorT, & ! Human activity
         TraffDorNorT, & ! Traffic
         AHDorNorT ! Anthropogenic heat
 
      ! initialisation
      qf_traff = 0
      qf_sahp_heating = 0
      qf_sahp_cooling = 0
 
      ! Transfer HDD/CDD values to local explict variables
      hdd_daily = hdd_id(7)
      cdd_daily = hdd_id(8)
      ! NB: temporarily use 5-day running mean to test the performance
      tair_avg_daily = hdd_id(10)
 
      ! Tair_avg_daily= HDD_id_use(3) ! this is daily
 
      ! use average of both daytime and nighttime population values
      ! TS 20 Sep 2019: moved from  here to simplify the interface
      IF (popdensdaytime(1) >= 0 .AND. popdensnighttime >= 0) numcapita(1) = (popdensdaytime(1) + popdensnighttime)/2 !If both, use average
      IF (popdensdaytime(2) >= 0 .AND. popdensnighttime >= 0) numcapita(2) = (popdensdaytime(2) + popdensnighttime)/2 !If both, use average
 
      !NumCapita = (PopDensDaytime + PopDensNighttime)/2
 
      !-----------------------------------------------------------------------
      ! Account for Daylight saving
      ih = it - dls
      IF (ih < 0) ih = 23
 
      ! Set weekday/weekend counter
      iu = 1 !Set to 1=weekday
      IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) iu = 2 !Set to 2=weekend
 
      ! Calculate energy emissions and CO2 from human metabolism -------------
      ! Pop dens (cap ha-1 -> cap m-2) x activity level (W cap-1)
      ! PopDorNorT   = PopProf_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)      ! 1=night, 2=day, 1-2=transition
      ! ActDorNorT   = HumActivity_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)  ! 1=night, 2=day, 1-2=transition
      ! TraffDorNorT = TraffProf_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)    ! normalise so the AVERAGE of the multipliers is equal to 1
      ! AHDorNorT    = AHProf_tstep((NSH*(ih+1-1)+imin*NSH/60+1),iu)       ! normalise so the AVERAGE of the multipliers is equal to 1
      ! PRINT*, ''
      ! PRINT*, 'AHDorNorT old:',AHDorNorT
      popdornort = get_prof_spectime_inst(ih, imin, 0, popprof_24hr(:, iu)) ! 1=night, 2=day, 1-2=transition
      actdornort = get_prof_spectime_inst(ih, imin, 0, humactivity_24hr(:, iu)) ! 1=night, 2=day, 1-2=transition
      traffdornort = get_prof_spectime_mean(ih, imin, 0, traffprof_24hr(:, iu)) ! normalise so the AVERAGE of the multipliers is equal to 1
      ahdornort = get_prof_spectime_mean(ih, imin, 0, ahprof_24hr(:, iu)) ! normalise so the AVERAGE of the multipliers is equal to 1
 
      ! Diurnal profile times population density [cap ha-1]
      dp_x_rhopop = ahdornort*numcapita(iu)
 
      qf_metab = (popdensnighttime*minqfmetab*((2 - actdornort) + (2 - popdornort))/2 + &
                  popdensdaytime(iu)*maxqfmetab*((actdornort - 1) + (popdornort - 1))/2)/10000 !W m-2
 
      fc_metab = (popdensnighttime*minfcmetab*((2 - actdornort) + (2 - popdornort))/2 + &
                  popdensdaytime(iu)*maxfcmetab*((actdornort - 1) + (popdornort - 1))/2)/10000 !umol m-2 s-1
 
      !----------------------------------------------------------------------------------------------------
      !Different methods to calculate the anthopogenic heat and CO2 emissions.
      !1-3: CO2 emission is calculated from QF which can be calculated with three methods
      !41-43: CO2 emission is calculated with local information. QF methods are used for housing and human metabolism
 
      IF (emissionsmethod == 1 .OR. emissionsmethod == 4 .OR. &
          emissionsmethod == 11 .OR. emissionsmethod == 14 .OR. &
          emissionsmethod == 21 .OR. emissionsmethod == 24 .OR. &
          emissionsmethod == 31 .OR. emissionsmethod == 34 .OR. &
          emissionsmethod == 41 .OR. emissionsmethod == 44) THEN ! (formerly SAHP_1 subroutine)
         ! Loridan et al. (2011) JAMC Eq 13: linear relation with air temperature
         ! Weekday/weekend differences due to profile only
         ! Now scales with population density
 
         IF (temp_c < baset_heating(iu)) THEN
            ! QF_SAHP = (AH_MIN(iu) + AH_SLOPE_Heating(iu)*(BaseT_Heating(iu) - Temp_C))*AHDorNorT
            qf_sahp_heating = (ah_slope_heating(iu)*(baset_heating(iu) - temp_c))*ahdornort
         ELSE
            qf_sahp_heating = 0
            !    QF_SAHP = AH_MIN(iu)*AHDorNorT
         END IF
 
         ! Need to be checked later, not recommended to use
         qf_sahp_base = ah_min(iu)*ahdornort ! Temperature-independent contribution
         qf_sahp_cooling = 0 ! No AC contribution with this method
 
      ELSEIF (emissionsmethod == 2 .OR. emissionsmethod == 5 .OR. &
              emissionsmethod == 12 .OR. emissionsmethod == 15 .OR. &
              emissionsmethod == 22 .OR. emissionsmethod == 25 .OR. &
              emissionsmethod == 32 .OR. emissionsmethod == 35 .OR. &
              emissionsmethod == 42 .OR. emissionsmethod == 45) THEN ! (formerly SAHP_2 subroutine)
         ! Jarvi et al. (2011) JH Eq 3 using HDD and CDD
         ! Weekday/weekend differences due to profile and coefficients QF_a,b,c
         ! Scales with population density
         ! 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!
         qf_sahp_base = (qf_a(iu))*dp_x_rhopop ! Temperature-independent contribution from buildings, traffic and human metabolism
         qf_sahp_heating = (qf_c(iu)*hdd_daily)*dp_x_rhopop ! Heating contribution
         qf_sahp_cooling = (qf_b(iu)*cdd_daily)*dp_x_rhopop ! Cooling (AC) contribution
 
      ELSEIF (emissionsmethod == 3 .OR. emissionsmethod == 6 .OR. &
              emissionsmethod == 13 .OR. emissionsmethod == 16 .OR. &
              emissionsmethod == 23 .OR. emissionsmethod == 26 .OR. &
              emissionsmethod == 33 .OR. emissionsmethod == 36 .OR. &
              emissionsmethod == 43 .OR. emissionsmethod == 46) THEN
         ! Updated Loridan et al. (2011) method using daily (not instantaneous) air temperature (HDD_id_use(3))
         ! Linear relation with air temperature
         ! Weekday/weekend differences due to profile only
         ! Scales with population density
 
         ! Need to be checked later, not recommended to use
         qf_sahp_base = ah_min(iu)*ahdornort ! Temperature-independent contribution
 
         IF (tair_avg_daily < baset_heating(iu)) THEN ! Heating
            qf_sahp_heating = (ah_slope_heating(iu)*(baset_heating(iu) - tair_avg_daily))*ahdornort ! Heating contribution
            qf_sahp_cooling = 0
 
         ELSEIF (tair_avg_daily > baset_cooling(iu)) THEN ! Air-conditioning
            qf_sahp_heating = 0
            qf_sahp_cooling = (ah_slope_cooling(iu)*(tair_avg_daily - baset_cooling(iu)))*ahdornort ! AC contribution
 
         ELSE
            qf_sahp_heating = 0
            qf_sahp_cooling = 0
         END IF
 
      END IF
 
      ! total QF as the sum of components
      qf_sahp = qf_sahp_base + qf_sahp_heating + qf_sahp_cooling
 
      IF (emissionsmethod >= 1 .AND. emissionsmethod <= 3 .OR. &
          emissionsmethod >= 11 .AND. emissionsmethod <= 13 .OR. &
          emissionsmethod >= 21 .AND. emissionsmethod <= 23 .OR. &
          emissionsmethod >= 31 .AND. emissionsmethod <= 33 .OR. &
          emissionsmethod >= 41 .AND. emissionsmethod <= 43) THEN
 
         ! Calculate QF from buildings.
         ! First remove (if possibe) human metabolism from the total value given by SAHP.
         IF ((qf_sahp_base - qf_metab) > 0) THEN
            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
            !relative to traffic as metabolism is separately calculated
         ELSE
            CALL errorhint(69, 'QF metab exceeds base QF.', qf_metab, qf_sahp_base)
 
            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
         END IF
 
         ! Consider the various components of QF_build to calculate Fc_build
         ! Assume all A/C electric, so QF_SAHP_ac is not associated with any local CO2 emissions
         ! HDD part is building energy use, split between electricity (no local emissions CO2) and combustion (CO2) heating...
         fc_build = qf_sahp_heating*frfossilfuel_heat*ef_umolco2perj
 
         ! ... and there is also a temperature-independent contribution from building energy use.
         IF ((qf_sahp_base - qf_metab) > 0) THEN
            fc_build = fc_build + qf_sahp_base*qf0_beu(iu)*frfossilfuel_nonheat*ef_umolco2perj
         END IF
 
         ! Remainder of temperature-independent part must be traffic
         !QF_traff = (QF_SAHP_base-QF_metab)*(1.0-QF0_BEU(iu))
         qf_traff = qf_sahp_base*(1.0 - qf0_beu(iu)) - qf_metab
 
         ! Divide QF by energy emission factor and multiply by CO2 factor
         fc_traff = qf_traff/enef_v_jkm*fcef_v_kgkm(iu)*1e3*1e6/44
 
         !CO2 point source [kg C day-1] * 1e3 g kg-1 * 1e6 umol mol-1 /(12 g mol-1 * m-2)
         IF (co2pointsource > 0) THEN
            fc_point = co2pointsource*1e3*1e6/(12*60*60*24*surfacearea)
         ELSE
            fc_point = 0
         END IF
 
         ! Sum components to give anthropogenic CO2 flux [umol m-2 s-1]
         fc_anthro = fc_metab + fc_traff + fc_build + fc_point
 
      ELSEIF (emissionsmethod >= 4 .AND. emissionsmethod <= 6 .OR. &
              emissionsmethod >= 14 .AND. emissionsmethod <= 16 .OR. &
              emissionsmethod >= 24 .AND. emissionsmethod <= 26 .OR. &
              emissionsmethod >= 34 .AND. emissionsmethod <= 36 .OR. &
              emissionsmethod >= 44 .AND. emissionsmethod <= 46) THEN
         ! Calculate QF and Fc using building energy use and transport statistics
 
         ! Calculate energy released from traffic
         IF (trafficunits == 1) THEN ! [veh km m-2 day-1]
            ! Calculate using mean traffic rate [veh km m-2 day-1] * emission factor [J km-1]
            qf_traff = trafficrate(iu)/(60*60*24)*enef_v_jkm*traffdornort
            ! 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
            fc_traff = trafficrate(iu)/(60*60*24)*fcef_v_kgkm(iu)*1e3*1e6/44*traffdornort
 
         ELSEIF (trafficunits == 2) THEN ! [veh km cap-1 day-1]
            dp_x_rhopop_traff = traffdornort*numcapita(iu)/10000 ! [cap m-2]
            ! Use mean traffic rate [veh km cap-1 day-1] * emission factor [J km-1]
            qf_traff = trafficrate(iu)/(60*60*24)*enef_v_jkm*dp_x_rhopop_traff
            ! 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
            fc_traff = trafficrate(iu)/(60*60*24)*fcef_v_kgkm(iu)*1e3*1e6/44*dp_x_rhopop_traff
 
         ELSE ! If TrafficUnits doesn't match possible units
            CALL errorhint(75, 'Check TrafficUnits', trafficunits, -999d1)
 
         END IF
 
         ! Energy released from buildings only
         ! Should buildings have their own profile? Now using population profile
         qf_build = ((qf_sahp_base*qf0_beu(iu) + qf_sahp_heating + qf_sahp_cooling)/dp_x_rhopop)* &
                    (popdensnighttime*(2 - popdornort) + popdensdaytime(iu)*(popdornort - 1))
 
         ! Consider the various components of QF_build to calculate Fc_build
         fc_build = qf_sahp_heating*frfossilfuel_heat*ef_umolco2perj
         ! ... and there is also a temperature-independent contribution from building energy use
         fc_build = fc_build + qf_sahp_base*qf0_beu(iu)*frfossilfuel_nonheat*ef_umolco2perj
 
         !CO2 point source [kg C day-1] * 1e3 g kg-1 * 1e6 umol mol-1 /(12 g mol-1 * m-2)
         IF (co2pointsource > 0) THEN
            fc_point = co2pointsource*1e3*1e6/(12*60*60*24*surfacearea)
         ELSE
            fc_point = 0
         END IF
 
         ! Add other QF components to QF_SAHP_base (because if AnthropEmissionsMethod = 4, QF calculated above is building part only)
         qf_sahp_base = qf_sahp_base + qf_traff + qf_metab
 
         ! Sum components to give anthropogenic heat flux [W m-2]
         qf_sahp = qf_metab + qf_traff + qf_build
 
         ! Sum components to give anthropogenic CO2 flux [umol m-2 s-1]
         fc_anthro = fc_metab + fc_traff + fc_build + fc_point
 
      END IF
 
      RETURN
   END SUBROUTINE anthropogenicemissions
 
END MODULE anemsn_module
!========================================================================================