Skip to content

Commit

Permalink
Refactor temperature variable names in suews_phys_anemsn.f95
Browse files Browse the repository at this point in the history
  • Loading branch information
sunt05 committed Dec 23, 2023
1 parent fee1955 commit 2d6b87f
Show file tree
Hide file tree
Showing 5 changed files with 127 additions and 277 deletions.
24 changes: 17 additions & 7 deletions src/suews/src/suews_ctrl_driver.f95
Original file line number Diff line number Diff line change
Expand Up @@ -472,9 +472,8 @@ SUBROUTINE SUEWS_cal_Main_DTS( &
IF (Diagnose == 1) WRITE (*, *) 'Calling NARP_cal_SunPosition...'
! print *, 'timer: azimuth, zenith_deg', timer%azimuth, timer%zenith_deg
CALL NARP_cal_SunPosition_DTS( &
timer, & !input:
siteInfo, &
azimuth_deg, zenith_deg) !output:
timer, config, forcing, siteInfo, & !input
solarState)

!=================Calculation of density and other water related parameters=================
IF (Diagnose == 1) WRITE (*, *) 'Calling LUMPS_cal_AtmMoist...'
Expand Down Expand Up @@ -557,6 +556,7 @@ SUBROUTINE SUEWS_cal_Main_DTS( &
CALL SUEWS_cal_AnthropogenicEmission_DTS( &
timer, config, forcing, siteInfo, & ! input
anthroEmisState, &
atmState, &
heatState)

! ========================================================================
Expand Down Expand Up @@ -2574,20 +2574,25 @@ END SUBROUTINE SUEWS_cal_AnthropogenicEmission
SUBROUTINE SUEWS_cal_AnthropogenicEmission_DTS( &
timer, config, forcing, siteInfo, & ! input
anthroEmisState, &
atmState, &
heatState)
! QF, &
! QF_SAHP, &
! Fc_anthro, Fc_build, Fc_metab, Fc_point, Fc_traff) ! output:

USE SUEWS_DEF_DTS, ONLY: SUEWS_SITE, SUEWS_TIMER, SUEWS_CONFIG, SUEWS_FORCING, &
anthroEmis_STATE
anthroEmis_STATE, atm_state

IMPLICIT NONE
TYPE(SUEWS_TIMER), INTENT(IN) :: timer
TYPE(SUEWS_CONFIG), INTENT(IN) :: config
TYPE(SUEWS_FORCING), INTENT(IN) :: forcing
TYPE(SUEWS_SITE), INTENT(IN) :: siteInfo

TYPE(anthroEmis_STATE), INTENT(INout) :: anthroEmisState
TYPE(heat_STATE), INTENT(INout) :: heatState
TYPE(atm_state), INTENT(INout) :: atmState

! INTEGER, INTENT(in)::Diagnose
INTEGER :: EmissionsMethod !0 - Use values in met forcing file, or default QF;1 - Method according to Loridan et al. (2011) : SAHP; 2 - Method according to Jarvi et al. (2011) : SAHP_2

Expand Down Expand Up @@ -2628,9 +2633,9 @@ SUBROUTINE SUEWS_cal_AnthropogenicEmission_DTS( &

INTEGER, PARAMETER :: notUsedI = -999
REAL(KIND(1D0)), PARAMETER :: notUsed = -999
REAL(KIND(1D0)) :: Temp_local ! local ambient air temperature [degC]


TYPE(anthroEmis_STATE), INTENT(INout) :: anthroEmisState
TYPE(heat_STATE), INTENT(INout) :: heatState
ASSOCIATE ( &
dayofWeek_id => timer%dayofWeek_id, &
DLS => timer%DLS, &
Expand All @@ -2639,6 +2644,7 @@ SUBROUTINE SUEWS_cal_AnthropogenicEmission_DTS( &

ASSOCIATE ( &
EmissionsMethod => config%EmissionsMethod, &
localClimateMethod => config%localClimateMethod, &
EF_umolCO2perJ => ahemisPrm%EF_umolCO2perJ, &
EnEF_v_Jkm => ahemisPrm%EnEF_v_Jkm, &
FcEF_v_kgkm => ahemisPrm%FcEF_v_kgkm, &
Expand All @@ -2661,6 +2667,7 @@ SUBROUTINE SUEWS_cal_AnthropogenicEmission_DTS( &
Fc_traff => anthroEmisState%Fc_traff, &
QF => heatState%QF, &
QF_SAHP => heatState%QF_SAHP, &
T2_c => atmstate%t2_C, &
Temp_C => forcing%Temp_C, &
QF_obs => forcing%QF_obs &
)
Expand Down Expand Up @@ -2703,14 +2710,17 @@ SUBROUTINE SUEWS_cal_AnthropogenicEmission_DTS( &
IF (EmissionsMethod == 0) THEN ! use observed qf
qf = QF_obs
ELSEIF ((EmissionsMethod > 0 .AND. EmissionsMethod <= 6) .OR. EmissionsMethod >= 11) THEN
! choose temperature for anthropogenic heat flux calculation
Temp_local = MERGE(T2_c, Temp_C, localClimateMethod == 1)

CALL 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, &
Temp_local, HDD_id, Qf_A, Qf_B, Qf_C, &
AH_MIN, AH_SLOPE_Heating, AH_SLOPE_Cooling, &
BaseT_Heating, BaseT_Cooling, &
TrafficRate, &
Expand Down
8 changes: 4 additions & 4 deletions src/suews/src/suews_phys_anemsn.f95
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ SUBROUTINE AnthropogenicEmissions( &
FrFossilFuel_Heat, FrFossilFuel_NonHeat, &
MinFCMetab, MaxFCMetab, MinQFMetab, MaxQFMetab, &
PopDensDaytime, PopDensNighttime, &
Temp_C, HDD_id, Qf_A, Qf_B, Qf_C, &
Temp_local, HDD_id, Qf_A, Qf_B, Qf_C, &
AH_MIN, AH_SLOPE_Heating, AH_SLOPE_Cooling, &
BaseT_Heating, BaseT_Cooling, &
TrafficRate, &
Expand Down Expand Up @@ -86,7 +86,7 @@ SUBROUTINE AnthropogenicEmissions( &
MaxQFMetab, & !Maximum QF Metabolism
PopDensNighttime, & !Nighttime population density [ha-1] (i.e. residents)
SurfaceArea, & !Surface area [m-2]
Temp_C !Air temperature
Temp_local !Local ambient air temperature [degC]

REAL(KIND(1D0)), INTENT(out) :: &
QF_SAHP, &
Expand Down Expand Up @@ -179,9 +179,9 @@ SUBROUTINE AnthropogenicEmissions( &
! Weekday/weekend differences due to profile only
! Now scales with population density

IF (Temp_C < BaseT_Heating(iu)) THEN
IF (Temp_local < 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
QF_SAHP_heating = (AH_SLOPE_Heating(iu)*(BaseT_Heating(iu) - Temp_local))*AHDorNorT
ELSE
QF_SAHP_heating = 0
! QF_SAHP = AH_MIN(iu)*AHDorNorT
Expand Down
7 changes: 0 additions & 7 deletions src/suews/src/suews_phys_beers.f95
Original file line number Diff line number Diff line change
Expand Up @@ -493,13 +493,6 @@ SUBROUTINE BEERS_cal_main_DTS( &
Diagnose => config%Diagnose &
)

! iy = timer%iy
! id = timer%id

! avkdn = forcing%kdown
! Temp_C = forcing%Temp_C
! avrh = forcing%RH
! Press_hPa = forcing%Pres
ASSOCIATE ( &
emis_ground => pavedPrm%emis, &
emis_wall => bldgPrm%emis, &
Expand Down
200 changes: 106 additions & 94 deletions src/suews/src/suews_phys_narp.f95
Original file line number Diff line number Diff line change
Expand Up @@ -626,20 +626,25 @@ SUBROUTINE NARP_cal_SunPosition(year, idectime, UTC, &

END SUBROUTINE NARP_cal_SunPosition

SUBROUTINE NARP_cal_SunPosition_DTS(timer, &
siteInfo, &
sunazimuth, sunzenith)
SUBROUTINE NARP_cal_SunPosition_DTS( &
timer, config, forcing, siteInfo, & !input
solarState)
! sunazimuth, sunzenith)

USE SUEWS_DEF_DTS, ONLY: SUEWS_TIMER, SUEWS_SITE
USE SUEWS_DEF_DTS, ONLY: &
SUEWS_TIMER, SUEWS_SITE, SUEWS_CONFIG, SUEWS_FORCING, &
solar_State

IMPLICIT NONE

TYPE(SUEWS_TIMER), INTENT(in) :: timer
TYPE(SUEWS_SITE), INTENT(in) :: siteInfo
! REAL(KIND(1D0)), INTENT(in) :: dectime
TYPE(SUEWS_TIMER), INTENT(IN) :: timer
TYPE(SUEWS_CONFIG), INTENT(IN) :: config
TYPE(SUEWS_FORCING), INTENT(IN) :: forcing
TYPE(SUEWS_SITE), INTENT(IN) :: siteInfo

REAL(KIND(1D0)) :: year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude
REAL(KIND(1D0)), INTENT(out) :: sunazimuth, sunzenith
TYPE(solar_State), INTENT(INOUT) :: solarState

REAL(KIND(1D0)) :: year, idectime

REAL(KIND(1D0)) :: sec
INTEGER :: month, day, hour, min, seas, dayofyear, year_int
Expand All @@ -661,92 +666,99 @@ SUBROUTINE NARP_cal_SunPosition_DTS(timer, &
REAL(KIND(1D0)) :: topocentric_sun_positiondeclination
REAL(KIND(1D0)) :: topocentric_local_hour

year = REAL(timer%iy, KIND(1D0))
idectime = timer%dectime - timer%tstep/2/86400
UTC = siteInfo%timezone
locationlatitude = siteInfo%lat
locationlongitude = siteInfo%lon
locationaltitude = siteInfo%alt

! REAL(KIND(1D0)) :: sunazimuth,sunzenith

! This function compute the sun position (zenith and azimuth angle (in degrees) at the observer
! location) as a function of the observer local time and position.
!
! Input lat and lng should be in degrees, alt in meters.
!
! It is an implementation of the algorithm presented by Reda et Andreas in:
! Reda, I., Andreas, A. (2003) Solar position algorithm for solar
! radiation application. National Renewable Energy Laboratory (NREL)
! Technical report NREL/TP-560-34302.
! This document is available at www.osti.gov/bridge
! Code is translated from matlab code by Fredrik Lindberg ([email protected])
! Last modified: LJ 27 Jan 2016 - Tabs removed

! Convert to timevectors from dectime and year
CALL dectime_to_timevec(idectime, hour, min, sec)
dayofyear = FLOOR(idectime)
year_int = INT(year)
CALL day2month(dayofyear, month, day, seas, year_int, locationlatitude)

! 1. Calculate the Julian Day, and Century. Julian Ephemeris day, century
! and millenium are calculated using a mean delta_t of 33.184 seconds.
CALL julian_calculation(year, month, day, hour, min, sec, UTC, juliancentury, julianday, julianephemeris_century, &
julianephemeris_day, julianephemeris_millenium)

! 2. Calculate the Earth heliocentric longitude, latitude, and radius
! vector (L, B, and R)
CALL earth_heliocentric_position_calculation(julianephemeris_millenium, earth_heliocentric_positionlatitude,&
&earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)

! 3. Calculate the geocentric longitude and latitude
CALL sun_geocentric_position_calculation(earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude,&
& sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)

! 4. Calculate the nutation in longitude and obliquity (in degrees).
CALL nutation_calculation(julianephemeris_century, nutationlongitude, nutationobliquity)

! 5. Calculate the true obliquity of the ecliptic (in degrees).
CALL corr_obliquity_calculation(julianephemeris_millenium, nutationobliquity, corr_obliquity)

! 6. Calculate the aberration correction (in degrees)
CALL abberation_correction_calculation(earth_heliocentric_positionradius, aberration_correction)

! 7. Calculate the apparent sun longitude in degrees)
CALL apparent_sun_longitude_calculation(sun_geocentric_positionlongitude, nutationlongitude,&
& aberration_correction, apparent_sun_longitude)

! 8. Calculate the apparent sideral time at Greenwich (in degrees)
CALL apparent_stime_at_greenwich_calculation(julianday, juliancentury, nutationlongitude, &
&corr_obliquity, apparent_stime_at_greenwich)

! 9. Calculate the sun rigth ascension (in degrees)
CALL sun_rigth_ascension_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
&sun_rigth_ascension)

! 10. Calculate the geocentric sun declination (in degrees). Positive or
! negative if the sun is north or south of the celestial equator.
CALL sun_geocentric_declination_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
&sun_geocentric_declination)

! 11. Calculate the observer local hour angle (in degrees, westward from south).
CALL observer_local_hour_calculation(apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)

! 12. Calculate the topocentric sun position (rigth ascension, declination and
! rigth ascension parallax in degrees)
CALL topocentric_sun_position_calculate(topocentric_sun_positionrigth_ascension,&
&topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination, locationaltitude,&
&locationlatitude, observer_local_hour, sun_rigth_ascension, sun_geocentric_declination,&
&earth_heliocentric_positionradius)

! 13. Calculate the topocentric local hour angle (in degrees)
CALL topocentric_local_hour_calculate(observer_local_hour, topocentric_sun_positionrigth_ascension_parallax,&
& topocentric_local_hour)

! 14. Calculate the topocentric zenith and azimuth angle (in degrees)
CALL sun_topocentric_zenith_angle_calculate(locationlatitude, topocentric_sun_positiondeclination,&
& topocentric_local_hour, sunazimuth, sunzenith)

ASSOCIATE ( &
iy => timer%iy, &
dectime => timer%dectime, &
tstep => timer%tstep, &
UTC => siteInfo%timezone, &
locationlatitude => siteInfo%lat, &
locationlongitude => siteInfo%lon, &
locationaltitude => siteInfo%alt, &
azimuth_deg => solarState%azimuth_deg, &
zenith_deg => solarState%zenith_deg &
)
year = REAL(iy, KIND(1D0))
idectime = dectime - tstep/2/86400

! This function compute the sun position (zenith and azimuth angle (in degrees) at the observer
! location) as a function of the observer local time and position.
!
! Input lat and lng should be in degrees, alt in meters.
!
! It is an implementation of the algorithm presented by Reda et Andreas in:
! Reda, I., Andreas, A. (2003) Solar position algorithm for solar
! radiation application. National Renewable Energy Laboratory (NREL)
! Technical report NREL/TP-560-34302.
! This document is available at www.osti.gov/bridge
! Code is translated from matlab code by Fredrik Lindberg ([email protected])
! Last modified: LJ 27 Jan 2016 - Tabs removed

! Convert to timevectors from dectime and year
CALL dectime_to_timevec(idectime, hour, min, sec)
dayofyear = FLOOR(idectime)
year_int = INT(year)
CALL day2month(dayofyear, month, day, seas, year_int, locationlatitude)

! 1. Calculate the Julian Day, and Century. Julian Ephemeris day, century
! and millenium are calculated using a mean delta_t of 33.184 seconds.
CALL julian_calculation(year, month, day, hour, min, sec, UTC, juliancentury, julianday, julianephemeris_century, &
julianephemeris_day, julianephemeris_millenium)

! 2. Calculate the Earth heliocentric longitude, latitude, and radius
! vector (L, B, and R)
CALL earth_heliocentric_position_calculation(julianephemeris_millenium, earth_heliocentric_positionlatitude,&
&earth_heliocentric_positionlongitude, earth_heliocentric_positionradius)

! 3. Calculate the geocentric longitude and latitude
CALL sun_geocentric_position_calculation(earth_heliocentric_positionlongitude, earth_heliocentric_positionlatitude,&
& sun_geocentric_positionlatitude, sun_geocentric_positionlongitude)

! 4. Calculate the nutation in longitude and obliquity (in degrees).
CALL nutation_calculation(julianephemeris_century, nutationlongitude, nutationobliquity)

! 5. Calculate the true obliquity of the ecliptic (in degrees).
CALL corr_obliquity_calculation(julianephemeris_millenium, nutationobliquity, corr_obliquity)

! 6. Calculate the aberration correction (in degrees)
CALL abberation_correction_calculation(earth_heliocentric_positionradius, aberration_correction)

! 7. Calculate the apparent sun longitude in degrees)
CALL apparent_sun_longitude_calculation(sun_geocentric_positionlongitude, nutationlongitude,&
& aberration_correction, apparent_sun_longitude)

! 8. Calculate the apparent sideral time at Greenwich (in degrees)
CALL apparent_stime_at_greenwich_calculation(julianday, juliancentury, nutationlongitude, &
&corr_obliquity, apparent_stime_at_greenwich)

! 9. Calculate the sun rigth ascension (in degrees)
CALL sun_rigth_ascension_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
&sun_rigth_ascension)

! 10. Calculate the geocentric sun declination (in degrees). Positive or
! negative if the sun is north or south of the celestial equator.
CALL sun_geocentric_declination_calculation(apparent_sun_longitude, corr_obliquity, sun_geocentric_positionlatitude, &
&sun_geocentric_declination)

! 11. Calculate the observer local hour angle (in degrees, westward from south).
CALL observer_local_hour_calculation( &
apparent_stime_at_greenwich, locationlongitude, sun_rigth_ascension, observer_local_hour)

! 12. Calculate the topocentric sun position (rigth ascension, declination and
! rigth ascension parallax in degrees)
CALL topocentric_sun_position_calculate(topocentric_sun_positionrigth_ascension,&
&topocentric_sun_positionrigth_ascension_parallax, topocentric_sun_positiondeclination, locationaltitude,&
&locationlatitude, observer_local_hour, sun_rigth_ascension, sun_geocentric_declination,&
&earth_heliocentric_positionradius)

! 13. Calculate the topocentric local hour angle (in degrees)
CALL topocentric_local_hour_calculate(observer_local_hour, topocentric_sun_positionrigth_ascension_parallax,&
& topocentric_local_hour)

! 14. Calculate the topocentric zenith and azimuth angle (in degrees)
CALL sun_topocentric_zenith_angle_calculate(locationlatitude, topocentric_sun_positiondeclination,&
& topocentric_local_hour, azimuth_deg, zenith_deg)
END ASSOCIATE
END SUBROUTINE NARP_cal_SunPosition_DTS

!================================ Subfunction definitions ========================================================!
Expand Down
Loading

0 comments on commit 2d6b87f

Please sign in to comment.