SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines | Variables
beers_module Module Reference

Functions/Subroutines

subroutine beers_cal_main (iy, id, dectime, lamdap, lamdaf, avkdn, ldown, temp_c, avrh, press_hpa, tsurf, lat, lng, alt, timezone, zenith_deg, azimuth, alb_ground, alb_bldg, emis_ground, emis_wall, dataoutlinebeers)
 
subroutine cal_ci_latenight (iy, doy, ta_degc, rh_frac, radg, lat, p_kpa, cilatenight, dectime_sunrise, zen_sunrise, i0_sunrise)
 
subroutine kroof (radi, radd, radg, f_sh, altitude, svfr, svfveg, shadow, psi, alb_bldg, kdown)
 
subroutine cal_svfalfa (svfr, svfveg, svfalfa, tmp)
 
real(kind(1d0)) function cal_ratio_height2width (lamdap, lamdaf)
 
real(kind(1d0)) function hwtosvf_ground (hw)
 
real(kind(1d0)) function hwtosvf_roof (hw)
 
subroutine tsurfbeers (iy, ta, rh, radi, i0, dectime, snup, altitude, zen, timezone, lat, lng, alt, tg, tgwall, altmax)
 
subroutine shadowgroundkusaka (hw, azimuth, zen, shadowground, shadowwalls)
 
subroutine clearnessindex_2013b (zen, doy, ta_degc, rh_frac, radg, lat, p_kpa, i0, ci, kt, i0et, ciuncorr)
 
subroutine sun_distance (jday, d)
 
subroutine cylindric_wedge (zen, svfalfa, f_sh)
 
subroutine diffusefraction (radg, altitude, kt, ta, rh, radi, radd)
 
subroutine kwalls (svf, svfveg, shadow, f_sh, radi, radg, radd, azimuth, altitude, psi, t, alb_ground, alb_bldg, keast, knorth, ksouth, kwest)
 
subroutine kvikt_veg (svf, svfveg, vikttot, viktveg, viktwall)
 
real(kind(1d0)) function cal_vikt (svf_x, vikttot)
 
subroutine lwalls (svf, svfveg, svfaveg, ldown2d, lup2d, altitude, ta, tw, sbc, emis_wall, emis_sky, t, ci, azimuth, ldown, svfalfa, f_sh_in, least, lnorth, lsouth, lwest)
 
subroutine lvikt_veg (isvf, isvfveg, isvfaveg, vikttot, viktveg, viktsky, viktrefl, viktwall)
 
subroutine issign (ix, maxpos, isignm)
 
subroutine day2month (b, mb, md, seas, year, latitude)
 
subroutine month2day (mon, ne, k, b)
 
subroutine leapyearcalc (year_int, nrodays)
 
elemental integer function days_of_year (year_int)
 
subroutine day_of_week (date, month, year, dow)
 
subroutine dectime_to_timevec (dectime, hours, mins, secs)
 
subroutine daylen (doy, xlat, dayl, dec, sndn, snup)
 
subroutine suews_cal_dectime (id, it, imin, isec, dectime)
 
subroutine suews_cal_tstep (tstep, nsh, nsh_real, tstep_real)
 
subroutine suews_cal_weekday (iy, id, lat, dayofweek_id)
 
subroutine suews_cal_dls (id, startdls, enddls, dls)
 

Variables

real(kind(1d0)), parameter pi = ATAN(1.)*4
 
real(kind(1d0)), parameter deg2rad = pi/180
 
real(kind(1d0)), parameter rad2deg = 1/DEG2RAD
 

Function/Subroutine Documentation

◆ beers_cal_main()

subroutine beers_module::beers_cal_main ( integer, intent(in) iy,
integer, intent(in) id,
real(kind(1d0)), intent(in) dectime,
real(kind(1d0)), intent(in) lamdap,
real(kind(1d0)), intent(in) lamdaf,
real(kind(1d0)), intent(in) avkdn,
real(kind(1d0)), intent(in) ldown,
real(kind(1d0)), intent(in) temp_c,
real(kind(1d0)), intent(in) avrh,
real(kind(1d0)), intent(in) press_hpa,
real(kind(1d0)), intent(in) tsurf,
real(kind(1d0)), intent(in) lat,
real(kind(1d0)), intent(in) lng,
real(kind(1d0)), intent(in) alt,
real(kind(1d0)), intent(in) timezone,
real(kind(1d0)), intent(in) zenith_deg,
real(kind(1d0)), intent(in) azimuth,
real(kind(1d0)), intent(in) alb_ground,
real(kind(1d0)), intent(in) alb_bldg,
real(kind(1d0)), intent(in) emis_ground,
real(kind(1d0)), intent(in) emis_wall,
real(kind(1d0)), dimension(ncolumnsdataoutbeers - 5), intent(out) dataoutlinebeers )

Definition at line 37 of file suews_phys_beers.f95.

42
43 IMPLICIT NONE
44
45 INTEGER, INTENT(in) :: iy
46 INTEGER, INTENT(in) :: id
47 REAL(KIND(1D0)), INTENT(in) :: lamdaP ! plan area fraction
48 REAL(KIND(1D0)), INTENT(in) :: lamdaF ! frontal area fraction
49 !REAL(KIND(1d0)), INTENT(in)::tilt ! Tilt of building in degrees. Could be included FL
50 ! REAL(KIND(1d0)), intent(in) ::lai_id_dectr
51 ! REAL(KIND(1d0)), intent(in) ::LAImax_dectr
52 ! REAL(KIND(1D0)), intent(in)::TransMin ! Tranmissivity of K through decidious vegetation (leaf on)
53 REAL(KIND(1D0)), INTENT(in) :: Press_hPa
54 REAL(KIND(1D0)), INTENT(in) :: Temp_C
55 REAL(KIND(1D0)), INTENT(in) :: avrh
56 REAL(KIND(1D0)), INTENT(in) :: avkdn
57 REAL(KIND(1D0)), INTENT(in) :: ldown
58 REAL(KIND(1D0)), INTENT(in) :: Tsurf
59 ! REAL(KIND(1d0)), intent(in) :: kdiff !Actual inputs from metfile.
60 ! REAL(KIND(1d0)), intent(in) :: kdir !Actual inputs from metfile.
61 REAL(KIND(1D0)), INTENT(in) :: zenith_deg
62 REAL(KIND(1D0)), INTENT(in) :: azimuth
63 REAL(KIND(1D0)), INTENT(in) :: dectime
64 REAL(KIND(1D0)), INTENT(in) :: timezone, lat, lng, alt
65 REAL(KIND(1D0)), INTENT(in) :: alb_bldg
66 REAL(KIND(1D0)), INTENT(in) :: alb_ground
67 REAL(KIND(1D0)), INTENT(in) :: emis_wall
68 REAL(KIND(1D0)), INTENT(in) :: emis_ground
69
70 REAL(KIND(1D0)), PARAMETER :: absL = 0.97 ! Absorption coefficient of longwave radiation of a person
71 REAL(KIND(1D0)), PARAMETER :: absK = 0.7 ! Absorption coefficient of shortwave radiation of a person
72 REAL(KIND(1D0)), PARAMETER :: Fside = 0.22 ! Standing human shape factor
73 REAL(KIND(1D0)), PARAMETER :: Fup = 0.06 ! Standing human shape factor
74
75 ! integer, parameter :: ncolumnsDataOutSol = 34 ! Standing human shape factor
76 REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutBEERS - 5), INTENT(OUT) :: dataOutLineBEERS ! 26 columns of output at the moment
77
78 REAL(KIND(1D0)) :: t, psi
79 REAL(KIND(1D0)) :: altitude, zen !azimuth,zenith
80 REAL(KIND(1D0)) :: CI, c, I0, Kt, Tw, Tg
81 REAL(KIND(1D0)) :: Ta, RH, P, radG, radD, radI !,idectime,tdectime!dectime,
82 REAL(KIND(1D0)) :: I0et, CIuncorr !,lati
83 REAL(KIND(1D0)) :: SNDN, SNUP, DEC, DAYL !,timestepdec,YEAR
84 REAL(KIND(1D0)) :: msteg, emis_sky, ea
85 REAL(KIND(1D0)) :: shadowground, shadowwalls, shadowroof
86 ! REAL(KIND(1d0)),intent(in) ::lai_id
87 INTEGER :: DOY !,ith!onlyglobal,usevegdem,x,y,i, first, second,
88 REAL(KIND(1D0)) :: CIlatenight
89 REAL(KIND(1D0)) :: dectime_sunrise, zen_sunrise, I0_sunrise
90 ! REAL(KIND(1d0)) :: Fside ! fraction of a person seen from each cardinal point
91 ! REAL(KIND(1d0)) :: Fup ! fraction of a person seen from down and up
92 REAL(KIND(1D0)) :: HW ! building height to width ratio
93
94 REAL(KIND(1D0)) :: svfalfa !, sos
95 !REAL(KIND(1d0)) :: gvf !Ground View Factors (GVF)
96 REAL(KIND(1D0)) :: Tmrt, Sstr, F_sh
97 !REAL(KIND(1d0)) :: vegsh
98 REAL(KIND(1D0)) :: tmp, altmax
99 REAL(KIND(1D0)) :: svf_bldg_veg
100 REAL(KIND(1D0)) :: svf_ground, svf_roof
101 REAL(KIND(1D0)) :: svf_veg
102 REAL(KIND(1D0)) :: svf_aveg
103 REAL(KIND(1D0)) :: Kdown, Keast, Knorth, Ksouth, Kup2d, Kwest
104 REAL(KIND(1D0)) :: Ldown2d, Least, Lnorth, Lsouth, Lup2d, Lwest
105
106 ! Internal grids
107 !Search directions for Ground View Factors (GVF)
108 !REAL(KIND(1d0)), PARAMETER :: azimuthA(1:18) = [(j*(360.0/18.0), j=0, 17)]
109 ! temporary parameters and variables for testing
110 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
111 REAL(KIND(1D0)), PARAMETER :: SBC = 5.67051e-8
112
113 INTEGER, PARAMETER :: onlyglobal = 1 !! force to calculate direct and diffuse components, TS 13 Dec 2019 !TODO: should be input parameter FL
114 !INTEGER, PARAMETER:: usevegdem = 0 !! force not to use vegetation DEM based calculations, TS 13 Dec 2019
115 !INTEGER, PARAMETER:: row = 1 !! force to 1, TS 13 Dec 2019
116 !INTEGER, PARAMETER:: col = 1 !! force to 1, TS 13 Dec 2019
117 !INTEGER, PARAMETER:: Posture = 1 !! force to 1, TS 13 Dec 2019 ! Not used FL
118 ! INTEGER, PARAMETER:: SOLWEIG_ldown = 0 !! force to 0, TS 13 Dec 2019 !TODO: should be input parameter FL
119 ! this is just for testing
120 INTEGER, PARAMETER :: SOLWEIG_ldown = 1 !! force to 0, TS 13 Dec 2019 !TODO: should be input parameter FL
121 INTEGER, PARAMETER :: BEERS_tsurf = 1 !TODO: should be input parameter FL
122
123!!!!!! Begin program !!!!!!
124 ! internal grids
125 ! ALLOCATE (tmp(1, 1))
126 !ALLOCATE (Knight(1, 1))
127 !ALLOCATE (Tgmap0(1, 1))
128 ! ALLOCATE (svfbuveg(1, 1))
129 !allocate(Tgmap0(1,1))
130
131 ! initialise this as ONE
132 cilatenight = 1
133 ci = 1
134
135 ! These variables should change name in the future...
136 p = press_hpa
137 ta = temp_c
138 rh = avrh
139 radg = avkdn
140 doy = int(id)
141 ! radD = kdiff
142 ! radI = kdir
143
144 psi = 0.03 ! Tranmissivity of K through vegetation
145
146 ! alb_bldg = alb(2) ! taken from Bldg (FunctionalTypes)
147 ! alb_ground = alb(1) ! taken from Paved (FunctionalTypes)
148 ! emis_wall = emis(2) ! taken from Bldg (FunctionalTypes)
149 ! emis_ground = emis(1) ! taken from Paved (FunctionalTypes)
150
151 ! Building azimuth offset from cardinal points in degrees
152 t = 0 !tilt
153
154 hw = cal_ratio_height2width(lamdap, lamdaf)
155
156 ! parameterisation using NYC building data
157 ! TODO: #5 which SVF should be used here? in python code, SVF_roof is used.
158 ! Both are used. svfr for roof fluxes. Changed in code below. FL
159 svf_ground = hwtosvf_ground(hw) !TODO: Should change based on Oke equation???
160 svf_roof = hwtosvf_roof(hw) !TODO: Should change based on Oke equation???
161
162 svf_veg = 1 ! svfveg: SVF based on vegetation blocking the sky (1 = no vegetation)
163 svf_aveg = 1 ! view factor where vegetation is in view before buildings.
164
165 tmp = 1 - (svf_ground + svf_veg - 1)
166 IF (tmp <= 1.e-6) tmp = 1.e-6 ! avoiding log(0)
167 svfalfa = asin(exp(log(tmp)/2))
168
169 ! SVF combines for buildings and vegetation
170 svf_bldg_veg = (svf_ground - (1 - svf_veg)*(1 - psi))
171
172 ! Sun position related things
173 CALL daylen(doy, lat, dayl, dec, sndn, snup)
174 zen = zenith_deg*deg2rad
175 altitude = 90 - zenith_deg
176
177 !Determination of clear-sky emissivity from Prata (1996) !TODO: Is this calcualted in NARP?
178 ea = 6.107*10**((7.5*ta)/(237.3 + ta))*(rh/100) !Vapor pressure
179 msteg = 46.5*(ea/(ta + 273.15))
180 emis_sky = (1 - (1 + msteg)*exp(-((1.2 + 3.0*msteg)**0.5)))
181
182 !!! DAYTIME !!!
183 IF (altitude > 0.1) THEN
184
185 !Clearness Index on Earth's surface after Crawford and Dunchon (1999) with a correction
186 !factor for low sun elevations after Lindberg et al. (2008)
187 CALL clearnessindex_2013b(zen, doy, ta, rh/100., radg, lat, p/10., & !input
188 i0, ci, kt, i0et, ciuncorr) !output
189 ! IF (CI > 1) CI = 1
190
191 !Estimation of radD and radI if not measured after Reindl et al. (1990)
192 IF (onlyglobal == 1) THEN
193 CALL diffusefraction(radg, altitude, kt, ta, rh, radi, radd)
194 END IF
195
196 CALL shadowgroundkusaka(hw, azimuth, zen, shadowground, shadowwalls)
197 shadowroof = 1. ! TODO: should change with time of day etc. Could be parameterizised from e.g. Lindberg et al. 2015 SE
198
199 CALL cylindric_wedge(zen, svfalfa, f_sh)
200
201 !!! Calculation of shortwave daytime radiative fluxes !!!
202 CALL kroof(radi, radd, radg, f_sh, altitude, svf_roof, svf_veg, shadowroof, psi, alb_bldg, kdown)
203 !Kdown2d = radI*shadowroof*SIN(altitude*DEG2RAD) &
204 ! + radD*svfr &
205 ! ! TODO: #6 F_sh issue: used below is calculated as a variable but
206 ! + alb_bldg*(1 - svfr)*(radG*(1 - F_sh) + radD*F_sh)
207
208 kup2d = alb_ground*( &
209 radi*shadowground*sin(altitude*deg2rad) & ! gvf not defined TODO #2 FIXED
210 + radd*svf_bldg_veg &
211 + alb_bldg*(1 - svf_bldg_veg)*(radg*(1 - f_sh) + radd*f_sh))
212
213 ! TODO: #7 check consistency with python code
214 CALL kwalls( &
215 svf_ground, svf_veg, shadowground, f_sh, &
216 radi, radg, radd, azimuth, altitude, psi, t, alb_ground, alb_bldg, & ! input
217 keast, knorth, ksouth, kwest) ! output
218
219 IF (beers_tsurf == 1) THEN
220 CALL tsurfbeers(iy, ta, rh, radi, i0, dectime, snup, altitude, zen, timezone, lat, lng, alt, &
221 tg, tw, altmax)
222 ELSE
223 tg = tsurf - ta !TODO: Tg is the difference (added temperature between TA and Tg) in BEERS
224 ! Tg = Tsurf ! changed to absolute sense
225 tw = tg
226 END IF
227
228 !!!! Lup, daytime !!!!
229 lup2d = sbc*emis_ground*((shadowground*tg + ta + 273.15)**4)
230
231 !!!!!!! NIGHTTIME !!!!!!!!
232 ELSE
233 CALL cal_ci_latenight(iy, doy, ta, rh/100., radg, lat, p/10., & !input
234 cilatenight, dectime_sunrise, zen_sunrise, i0_sunrise) !output
235 ci = cilatenight
236 i0 = i0_sunrise
237 shadowground = 0
238 shadowwalls = 0
239 shadowroof = 0
240 ! Tg = Temp_C !TODO: #11 This need some thought. Use ESTM to improve?
241 ! Tw = Temp_C !TODO: This need some thought. Use ESTM to improve?
242 IF (beers_tsurf == 1) THEN
243 tw = 0
244 tg = 0
245 ELSE
246 tg = tsurf - ta !TODO: Tg is the difference (added temperature between Ta and Tg) in BEERS
247 tw = tg
248 END IF
249 radi = 0
250 radd = 0
251 f_sh = 1
252
253 !Nocturnal cloud fraction from Offerle et al. 2003
254 IF (dectime < (doy + 0.5) .AND. dectime > doy .AND. altitude < 1.0) THEN !TODO: THIS (STILL, 20201117) NEED SOME THOUGHT 20150211
255 !j=0
256 !do while (dectime<(DOY+SNUP/24))
257 !! call ConvertMetData(ith+j) ! read data at sunrise ??
258 ! j=j+1
259 !end do
260 !call NARP_cal_SunPosition(year,dectime,timezone,lat,lng,alt,azimuth,zenith_deg)!this is not good
261 !zen=zenith_deg*DEG2RAD
262 !call clearnessindex_2013b(zen,DOY,Temp_C,RH/100,avkdn,lat,Press_hPa,I0,CI,Kt,I0et,CIuncorr)
263 !!call ConvertMetData(ith) ! read data at current timestep again ??
264 !call NARP_cal_SunPosition(year,dectime,timezone,lat,lng,alt,azimuth,zenith_deg)!this is not good
265
266 ci = 1.0
267 ELSE
268 ! IF (SolweigCount == 1) THEN
269 ! CI = 1.0
270 ! ELSE
271 ! CI = CIlatenight
272 ! END IF
273 ci = cilatenight
274 END IF
275
276 radi = 0
277 radd = 0
278
279 !Nocturnal Kfluxes set to 0
280 kdown = 0.0
281 kwest = 0.0
282 kup2d = 0.0
283 keast = 0.0
284 ksouth = 0.0
285 knorth = 0.0
286
287 !!! Lup !!!
288 lup2d = sbc*emis_ground*((ta + tg + 273.15)**4)
289
290 END IF
291
292 !!! Ldown !!!
293 IF (solweig_ldown == 1) THEN ! Third
294 ldown2d = (svf_roof + svf_veg - 1)*emis_sky*sbc*((ta + 273.15)**4) &
295 + (2 - svf_veg - svf_aveg)*emis_wall*sbc*((ta + 273.15)**4) &
296 + (svf_aveg - svf_roof)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
297 + (2 - svf_roof - svf_veg)*(1 - emis_wall)*emis_sky*sbc*((ta + 273.15)**4)
298
299 IF (ci < 0.95) THEN !non-clear conditions
300 c = 1 - ci
301 ldown2d = ldown2d*(1 - c) + c*( &
302 (svf_roof + svf_veg - 1)*sbc*((ta + 273.15)**4) &
303 + (2 - svf_veg - svf_aveg)*emis_wall*sbc*((ta + 273.15)**4) &
304 + (svf_aveg - svf_roof)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
305 + (2 - svf_roof - svf_veg)*(1 - emis_wall)*sbc*((ta + 273.15)**4))
306 END IF
307
308 ELSE
309 ldown2d = (svf_roof + svf_veg - 1)*ldown &
310 + (2 - svf_veg - svf_aveg)*emis_wall*sbc*((ta + 273.15)**4) &
311 + (svf_aveg - svf_roof)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
312 + (2 - svf_roof - svf_veg)*(1 - emis_wall)*ldown
313 END IF
314
315 !!! Lside !!!
316 CALL lwalls(svf_ground, svf_veg, svf_aveg, &
317 ldown2d, lup2d, &
318 altitude, ta, tw, sbc, emis_wall, &
319 emis_sky, t, ci, azimuth, ldown, svfalfa, f_sh, &
320 least, lnorth, lsouth, lwest) ! output
321
322 !!! Calculation of radiant flux density and Tmrt (Mean radiant temperature) !!!
323 sstr = absk*(kdown*fup + kup2d*fup + knorth*fside + keast*fside + ksouth*fside + kwest*fside) &
324 + absl*(ldown2d*fup + lup2d*fup + lnorth*fside + least*fside + lsouth*fside + lwest*fside)
325 tmrt = sqrt(sqrt((sstr/(absl*sbc)))) - 273.15
326
327 dataoutlinebeers = [azimuth, altitude, radg, radi, radd, &
328 kdown, kup2d, ksouth, kwest, knorth, keast, &
329 ldown2d, lup2d, lsouth, lwest, lnorth, least, &
330 tmrt, i0, ci, shadowground, shadowwalls, svf_ground, svf_roof, svf_bldg_veg, &
331 emis_sky, &
332 ta, tg, tw]
333
334 ! DEALLOCATE (tmp)
335 ! DEALLOCATE (svfbuveg)
336
daylen
subroutine daylen(doy, xlat, dayl, dec, sndn, snup)
Definition suews_util_time.f95:227

References cal_ci_latenight(), cal_ratio_height2width(), clearnessindex_2013b(), cylindric_wedge(), daylen(), allocatearray::deg2rad, diffusefraction(), hwtosvf_ground(), hwtosvf_roof(), kroof(), kwalls(), lwalls(), shadowgroundkusaka(), and tsurfbeers().

Referenced by suews_driver::suews_cal_main(), and suews_driver::suews_cal_main_dts().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_ci_latenight()

subroutine beers_module::cal_ci_latenight ( integer, intent(in) iy,
integer, intent(in) doy,
real(kind(1d0)), intent(in) ta_degc,
real(kind(1d0)), intent(in) rh_frac,
real(kind(1d0)), intent(in) radg,
real(kind(1d0)), intent(in) lat,
real(kind(1d0)), intent(in) p_kpa,
real(kind(1d0)), intent(out) cilatenight,
real(kind(1d0)), intent(out) dectime_sunrise,
real(kind(1d0)), intent(out) zen_sunrise,
real(kind(1d0)), intent(out) i0_sunrise )

Definition at line 339 of file suews_phys_beers.f95.

341 ! subroutine to calculate nighttime clearness index
342 ! Offerle et al. (2003) used sunset value but here a value after sunrise is calculated
343 ! TODO: a value at sunset can be retained using module variables
344
345 IMPLICIT NONE
346 INTEGER, INTENT(in) :: iy
347 INTEGER, INTENT(in) :: DOY
348 ! INTEGER, intent(in) :: timezone
349 REAL(KIND(1D0)), INTENT(in) :: Ta_degC
350 REAL(KIND(1D0)), INTENT(in) :: RH_frac
351 REAL(KIND(1D0)), INTENT(in) :: P_kPa
352 REAL(KIND(1D0)), INTENT(in) :: radG
353 REAL(KIND(1D0)), INTENT(in) :: lat
354 REAL(KIND(1D0)), INTENT(out) :: CIlatenight
355 REAL(KIND(1D0)), INTENT(out) :: dectime_sunrise
356 REAL(KIND(1D0)), INTENT(out) :: zen_sunrise
357 REAL(KIND(1D0)), INTENT(out) :: I0_sunrise
358 REAL(KIND(1D0)) :: I0et
359 REAL(KIND(1D0)) :: CIuncorr
360 REAL(KIND(1D0)) :: Kt
361 REAL(KIND(1D0)) :: DAYL, DEC, SNDN, SNUP, azimuth
362
363 CALL daylen(doy, lat, dayl, dec, sndn, snup)
364 dectime_sunrise = doy + (snup + .5)/24.
365 CALL narp_cal_sunposition( &
366 REAL(iy, KIND(1D0)), & !input:
367 dectime_sunrise, & ! sun position at middle of timestep before
368 0.d0, lat, 0.d0, 100.d0, &
369 azimuth, zen_sunrise) !output:
370
371 zen_sunrise = zen_sunrise/180.*pi
372
373 CALL clearnessindex_2013b(zen_sunrise, doy, ta_degc, rh_frac, radg, lat, p_kpa, & !input
374 i0_sunrise, cilatenight, kt, i0et, ciuncorr) !output
375

References clearnessindex_2013b(), daylen(), narp_module::narp_cal_sunposition(), and pi.

Referenced by beers_cal_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cal_ratio_height2width()

real(kind(1d0)) function beers_module::cal_ratio_height2width ( real(kind(1d0)), intent(in) lamdap,
real(kind(1d0)), intent(in) lamdaf )

Definition at line 438 of file suews_phys_beers.f95.

439 IMPLICIT NONE
440 REAL(KIND(1D0)), PARAMETER :: a = 0.5598
441 REAL(KIND(1D0)), PARAMETER :: b = -0.2485
442 REAL(KIND(1D0)), PARAMETER :: c = 0.4112
443 REAL(KIND(1D0)), PARAMETER :: d = -0.02528
444
445 REAL(KIND(1D0)), INTENT(in) :: lamdaP ! plan area fraction
446 REAL(KIND(1D0)), INTENT(in) :: lamdaF ! frontal area fraction
447 REAL(KIND(1D0)) :: HW
448 REAL(KIND(1D0)) :: lamdaW !wall area fraction (wallarea / total area)
449
450 ! TODO: we need to use wall area directly to get this
451 lamdaw = 4*lamdaf ! assuming square shaped buildings
452
453 hw = (lamdaw*lamdap)/(2*lamdap*(1 - lamdap))
454

Referenced by beers_cal_main().

Here is the caller graph for this function:

◆ cal_svfalfa()

subroutine beers_module::cal_svfalfa ( real(kind(1d0)), intent(in) svfr,
real(kind(1d0)), intent(in) svfveg,
real(kind(1d0)), intent(out) svfalfa,
real(kind(1d0)), intent(out) tmp )

Definition at line 422 of file suews_phys_beers.f95.

423 IMPLICIT NONE
424
425 REAL(KIND(1D0)), INTENT(in) :: svfr ! svf of roof
426 REAL(KIND(1D0)), INTENT(in) :: svfveg ! svf of vegetation
427 REAL(KIND(1D0)), INTENT(out) :: svfalfa ! Building height angle from svfr
428 REAL(KIND(1D0)), INTENT(out) :: tmp
429
430 ! TODO: we need to use wall area directly to get this
431 !Building height angle from svfr
432 tmp = 1.-(svfr + svfveg - 1.)
433 IF (tmp <= 0) tmp = 0.000000001 ! avoiding log(0)
434 svfalfa = asin(exp(log(tmp)/2.))
435

◆ cal_vikt()

real(kind(1d0)) function beers_module::cal_vikt ( real(kind(1d0)), intent(in) svf_x,
real(kind(1d0)), intent(in) vikttot )

Definition at line 949 of file suews_phys_beers.f95.

950 IMPLICIT NONE
951 REAL(KIND(1D0)), INTENT(in) :: svf_x, vikttot
952 REAL(KIND(1D0)) :: vikt
953
954 vikt = (vikttot &
955 - (63.227*svf_x**6 - 161.51*svf_x**5 &
956 + 156.91*svf_x**4 - 70.424*svf_x**3 &
957 + 16.773*svf_x**2 - 0.4863*svf_x))/vikttot
958

Referenced by kvikt_veg().

Here is the caller graph for this function:

◆ clearnessindex_2013b()

subroutine beers_module::clearnessindex_2013b ( real(kind(1d0)), intent(in) zen,
integer, intent(in) doy,
real(kind(1d0)), intent(in) ta_degc,
real(kind(1d0)), intent(in) rh_frac,
real(kind(1d0)), intent(in) radg,
real(kind(1d0)), intent(in) lat,
real(kind(1d0)), intent(in) p_kpa,
real(kind(1d0)), intent(out) i0,
real(kind(1d0)), intent(out) ci,
real(kind(1d0)), intent(out) kt,
real(kind(1d0)), intent(out) i0et,
real(kind(1d0)), intent(out) ciuncorr )

Definition at line 585 of file suews_phys_beers.f95.

587 ! Clearness Index at the Earth's surface calculated from Crawford and Duchon 1999
588 IMPLICIT NONE
589
590 INTEGER, INTENT(in) :: DOY
591 REAL(KIND(1D0)), INTENT(in) :: zen
592 REAL(KIND(1D0)), INTENT(in) :: Ta_degC
593 REAL(KIND(1D0)), INTENT(in) :: RH_frac
594 REAL(KIND(1D0)), INTENT(in) :: P_kPa
595 REAL(KIND(1D0)), INTENT(in) :: radG
596 REAL(KIND(1D0)), INTENT(in) :: lat
597 REAL(KIND(1D0)), INTENT(out) :: I0
598 REAL(KIND(1D0)), INTENT(out) :: CI
599 REAL(KIND(1D0)), INTENT(out) :: Kt
600 REAL(KIND(1D0)), INTENT(out) :: I0et
601 REAL(KIND(1D0)), INTENT(out) :: CIuncorr
602 REAL(KIND(1D0)) :: iG, Itoa, p
603 REAL(KIND(1D0)), DIMENSION(4) :: G
604 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
605
606 ! Variable declarations
607 real*8 :: a2
608 !logical :: b !>
609 real*8 :: b2
610 real*8 :: corr
611 real*8 :: d
612 real*8 :: m
613 real*8 :: tar
614 real*8 :: td
615 real*8 :: trpg
616 real*8 :: tw
617 real*8 :: u
618
619 IF (p_kpa == -999) THEN
620 p = 1013 !Pressure in millibars
621 ELSE
622 p = p_kpa*10 !Convert from hPa to millibars
623 END IF
624 !Effective solar constant
625 itoa = 1370 !TODO: Check if this is defined somewhere else?
626 CALL sun_distance(doy, d)
627 m = 35.*cos(zen)*((1224.*(cos(zen)**2) + 1.)**(-1./2.)) !optical air mass at p=1013
628 trpg = 1.021 - 0.084*(m*(0.000949*p + 0.051))**0.5 !Transmission coefficient for Rayliegh scattering and permanent gases
629
630 ! empirical constant depending on latitude
631 g = 0
632 IF (lat < 10) THEN
633 g = [3.37, 2.85, 2.80, 2.64]
634 ELSE IF (lat >= 10 .AND. lat < 20) THEN
635 g = [2.99, 3.02, 2.70, 2.93]
636 ELSE IF (lat >= 20 .AND. lat < 30) THEN
637 g = [3.60, 3.00, 2.98, 2.93]
638 ELSE IF (lat >= 30 .AND. lat < 40) THEN
639 g = [3.04, 3.11, 2.92, 2.94]
640 ELSE IF (lat >= 40 .AND. lat < 50) THEN
641 g = [2.70, 2.95, 2.77, 2.71]
642 ELSE IF (lat >= 50 .AND. lat < 60) THEN
643 g = [2.52, 3.07, 2.67, 2.93]
644 ELSE IF (lat >= 60 .AND. lat < 70) THEN
645 g = [1.76, 2.69, 2.61, 2.61]
646 ELSE IF (lat >= 70 .AND. lat < 80) THEN
647 g = [1.60, 1.67, 2.24, 2.63]
648 ELSE IF (lat >= 80 .AND. lat < 90) THEN
649 g = [1.11, 1.44, 1.94, 2.02]
650 END IF
651 IF (doy > 335 .OR. doy <= 60) THEN
652 ig = g(1)
653 ELSE IF (doy > 60 .AND. doy <= 152) THEN
654 ig = g(2)
655 ELSE IF (doy > 152 .AND. doy <= 244) THEN
656 ig = g(3)
657 ELSE IF (doy > 244 .AND. doy <= 335) THEN
658 ig = g(4)
659 END IF
660 !dewpoint calculation
661 a2 = 17.27
662 b2 = 237.7
663 td = (b2*(((a2*ta_degc)/(b2 + ta_degc)) + log(rh_frac)))/(a2 - (((a2*ta_degc)/(b2 + ta_degc)) + log(rh_frac)))
664 td = (td*1.8) + 32 !Dewpoint (degF)
665 u = exp(0.1133 - log(ig + 1) + 0.0393*td) !Precipitable water
666 tw = 1 - 0.077*((u*m)**0.3) !Transmission coefficient for water vapor
667 tar = 0.935**m !Transmission coefficient for aerosols
668
669 i0 = itoa*cos(zen)*trpg*tw*d*tar
670
671!!! This needs to be checked in fortran environment!!!
672 !b=I0==abs(zen)>pi/2
673 !I0(b==1)=0
674 !clear b
675 !if (not(isreal(I0))) then
676 ! I0=0
677 !end if
678
679 corr = 0.1473*log(90 - (zen/pi*180)) + 0.3454 ! 20070329
680
681 ciuncorr = radg/i0
682 ci = ciuncorr + (1 - corr)
683 i0et = itoa*cos(zen)*d !extra terrestial solar radiation
684 kt = radg/i0et
685
686 IF (ci > 1) ci = 1
687

References allocatearray::a2, pi, and sun_distance().

Referenced by beers_cal_main(), and cal_ci_latenight().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cylindric_wedge()

subroutine beers_module::cylindric_wedge ( real(kind(1d0)), intent(in) zen,
real(kind(1d0)), intent(in) svfalfa,
real(kind(1d0)), intent(out) f_sh )

Definition at line 706 of file suews_phys_beers.f95.

707 ! Fraction of sunlit walls based on sun altitude and svf wieghted building angles
708
709 IMPLICIT NONE
710
711 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
712 REAL(KIND(1D0)), INTENT(in) :: zen
713 REAL(KIND(1D0)), INTENT(in) :: svfalfa
714 REAL(KIND(1D0)), INTENT(out) :: F_sh
715 REAL(KIND(1D0)) :: beta
716 REAL(KIND(1D0)) :: alfa, xa, ha, hkil, ba
717 REAL(KIND(1D0)) :: Ai, phi, qa, Za
718 REAL(KIND(1D0)) :: ukil, Ssurf
719
720 ! ALLOCATE (alfa(1, 1))
721 ! ALLOCATE (ba(1, 1))
722 ! ALLOCATE (ha(1, 1))
723 ! ALLOCATE (xa(1, 1))
724 ! ALLOCATE (qa(1, 1))
725 ! ALLOCATE (Za(1, 1))
726 ! ALLOCATE (phi(1, 1))
727 ! ALLOCATE (ukil(1, 1))
728 ! ALLOCATE (Ai(1, 1))
729 ! ALLOCATE (Ssurf(1, 1))
730 ! ALLOCATE (hkil(1, 1))
731
732 beta = zen
733 alfa = svfalfa
734
735 xa = 1.-2./(tan(alfa)*tan(beta))
736 ha = 2./(tan(alfa)*tan(beta))
737 ba = (1./tan(alfa))
738 hkil = 2.*ba*ha
739 qa = 0.0d0
740
741 IF (xa < 0) qa = tan(beta)/2
742
743 za = 0.0d0
744 phi = 0.0d0
745 ai = 0.0d0
746 ukil = 0.0d0
747 IF (xa < 0) THEN
748 za = (ba**2 - qa**2/4.)**0.5
749 phi = atan(za/qa)
750 ai = (sin(phi) - phi*cos(phi))/(1 - cos(phi))
751 ukil = 2*ba*xa*ai
752 END IF
753
754 ssurf = hkil + ukil
755
756 f_sh = (2*pi*ba - ssurf)/(2*pi*ba) !Xa
757
758 IF (f_sh < 0) f_sh = 0.0
759 IF (f_sh > 0.5) f_sh = 0.5
760
761 ! DEALLOCATE (alfa)
762 ! DEALLOCATE (ba)
763 ! DEALLOCATE (ha)
764 ! DEALLOCATE (xa)
765 ! DEALLOCATE (qa)
766 ! DEALLOCATE (Za)
767 ! DEALLOCATE (phi)
768 ! DEALLOCATE (ukil)
769 ! DEALLOCATE (Ai)
770 ! DEALLOCATE (Ssurf)
771 ! DEALLOCATE (hkil)
772

References pi.

Referenced by beers_cal_main().

Here is the caller graph for this function:

◆ day2month()

subroutine beers_module::day2month ( integer, intent(in) b,
integer, intent(out) mb,
integer, intent(out) md,
integer, intent(out) seas,
integer, intent(in) year,
real(kind(1d0)) latitude )

Definition at line 1266 of file suews_phys_beers.f95.

1267 IMPLICIT NONE
1268 INTEGER, INTENT(in) :: b !b=doy --IN
1269 INTEGER, INTENT(out) :: mb !month=mb --OUT
1270 INTEGER, INTENT(out) :: md !date=md --OUT
1271 INTEGER, INTENT(out) :: seas
1272 INTEGER, INTENT(in) :: year
1273 INTEGER :: t1, t2, t3
1274 INTEGER :: k ! k- accounts for leap year
1275
1276 REAL(KIND(1D0)) :: latitude
1277
1278 ! initialisation
1279 mb = 1
1280
1281 !Corrected and calculation of date added LJ (Jun 2010)
1282
1283 t1 = 4
1284 t2 = 100
1285 t3 = 400
1286
1287 IF ((modulo(year, t1) == 0) .AND. (modulo(year, t2) /= 0) .OR. (modulo(year, t3) == 0)) THEN
1288 k = 1
1289 ELSE
1290 k = 0
1291 END IF
1292
1293 IF (b <= 31) THEN !January
1294 mb = 1
1295 md = b
1296 ELSEIF (b > 31 .AND. b <= 59 + k) THEN
1297 mb = 2
1298 md = b - 31
1299 ELSEIF (b > 59 + k .AND. b <= 90 + k) THEN
1300 mb = 3
1301 md = b - (59 + k)
1302 ELSEIF (b > 90 + k .AND. b <= 120 + k) THEN
1303 mb = 4
1304 md = b - (90 + k)
1305 ELSEIF (b > 120 + k .AND. b <= 151 + k) THEN
1306 mb = 5
1307 md = b - (120 + k)
1308 ELSEIF (b > 151 + k .AND. b <= 181 + k) THEN
1309 mb = 6
1310 md = b - (151 + k)
1311 ELSEIF (b > 181 + k .AND. b <= 212 + k) THEN
1312 mb = 7
1313 md = b - (181 + k)
1314 ELSEIF (b > 212 + k .AND. b <= 243 + k) THEN
1315 mb = 8
1316 md = b - (212 + k)
1317 ELSEIF (b > 243 + k .AND. b <= 273 + k) THEN
1318 mb = 9
1319 md = b - (243 + k)
1320 ELSEIF (b > 273 + k .AND. b <= 304 + k) THEN
1321 mb = 10
1322 md = b - (273 + k)
1323 ELSEIF (b > 304 + k .AND. b <= 334 + k) THEN
1324 mb = 11
1325 md = b - (304 + k)
1326 ELSEIF (b > 334 + k) THEN
1327 mb = 12
1328 md = b - (334 + k)
1329 END IF
1330
1331 !
1332 IF (latitude > 0) THEN ! Northern Hemisphere
1333 IF (mb > 3 .AND. mb < 10) THEN !Summer is from Apr to Sep
1334 seas = 1
1335 ELSE
1336 seas = 2 !Winter rest of the months
1337 END IF
1338 ELSE ! southern hemisphere
1339 IF (mb < 4 .OR. mb > 9) THEN !Summer is from Oct to Mar
1340 seas = 1
1341 ELSE
1342 seas = 2 !Winter rest of the months
1343 END IF
1344 END IF
1345 RETURN

Referenced by suews_cal_weekday().

Here is the caller graph for this function:

◆ day_of_week()

subroutine beers_module::day_of_week ( integer date,
integer month,
integer year,
integer dow )

Definition at line 1416 of file suews_phys_beers.f95.

1417 ! Calculate weekday from year, month and day information.
1418 ! DOW: Sunday=1,...Saturday=7
1419 ! YEAR fixed to integer, LJ March 2015
1420
1421 IMPLICIT NONE
1422
1423 INTEGER DATE, MONTH, DAY, YR, MN, N1, N2, DOW, YEAR
1424
1425 yr = year
1426 mn = month
1427
1428 !C
1429 !C IF JANUARY OR FEBRUARY, ADJUST MONTH AND YEAR
1430 !C
1431 IF (mn > 2) GO TO 10
1432 mn = mn + 12
1433 yr = yr - 1
143410 n1 = (26*(mn + 1))/10
1435 n2 = (125*yr)/100
1436 day = (date + n1 + n2 - (yr/100) + (yr/400) - 1)
1437 dow = mod(day, 7) + 1
1438
1439 RETURN

Referenced by suews_cal_weekday().

Here is the caller graph for this function:

◆ daylen()

subroutine beers_module::daylen ( integer doy,
real(kind(1d0)), intent(in) xlat,
real(kind(1d0)), intent(out) dayl,
real(kind(1d0)), intent(out) dec,
real(kind(1d0)), intent(out) sndn,
real(kind(1d0)), intent(out) snup )

Definition at line 1469 of file suews_phys_beers.f95.

1470 !=======================================================================
1471 ! DAYLEN, Real Function, N.B. Pickering, 09/23/1993
1472 ! Computes solar day length (Spitters, 1986).
1473 !=======================================================================
1474 !-----------------------------------------------------------------------
1475 IMPLICIT NONE
1476 INTEGER :: DOY
1477 REAL(KIND(1D0)), INTENT(IN) :: XLAT
1478 REAL(KIND(1D0)), INTENT(OUT) :: DEC, DAYL, SNDN, SNUP
1479 REAL(KIND(1D0)) :: SOC
1480 REAL(KIND(1D0)), PARAMETER :: RAD = pi/180.0
1481
1482 !-----------------------------------------------------------------------
1483 ! Calculation of declination of sun (Eqn. 16). Amplitude= +/-23.45
1484 ! deg. Minimum = DOY 355 (DEC 21), maximum = DOY 172.5 (JUN 21/22).
1485 dec = -23.45*cos(2.0*pi*(doy + 10.0)/365.0)
1486
1487 ! Sun angles. SOC limited for latitudes above polar circles.
1488 soc = tan(rad*dec)*tan(rad*xlat)
1489 soc = min(max(soc, -1.0), 1.0)
1490
1491 ! Calculate daylength, sunrise and sunset (Eqn. 17)
1492 dayl = 12.0 + 24.0*asin(soc)/pi
1493 snup = 12.0 - dayl/2.0
1494 sndn = 12.0 + dayl/2.0
1495

References pi.

Referenced by beers_cal_main(), and cal_ci_latenight().

Here is the caller graph for this function:

◆ days_of_year()

elemental integer function beers_module::days_of_year ( integer, intent(in) year_int)

Definition at line 1399 of file suews_phys_beers.f95.

1400 IMPLICIT NONE
1401 INTEGER, INTENT(in) :: year_int
1402 INTEGER :: nDays
1403
1404 IF (mod(year_int, 100) /= 0 .AND. mod(year_int, 4) == 0) THEN
1405 ndays = 366
1406 ELSEIF (mod(year_int, 400) == 0) THEN
1407 ndays = 366
1408 ELSE
1409 ndays = 365
1410 END IF
1411

References allocatearray::ndays.

◆ dectime_to_timevec()

subroutine beers_module::dectime_to_timevec ( real(kind(1d0)) dectime,
integer hours,
integer mins,
real(kind(1d0)) secs )

Definition at line 1445 of file suews_phys_beers.f95.

1446 !This subroutine converts dectime to individual
1447 !hours, minutes and seconds
1448 INTEGER :: HOURS, MINS, doy
1449 REAL(KIND(1D0)) :: dectime, SECS, DH, DM, DS
1450 !INTEGER :: year
1451
1452 doy = floor(dectime)
1453
1454 dh = dectime - doy !Decimal hours
1455 hours = int(24*dh)
1456
1457 dm = 24*dh - hours !Decimal minutes
1458 mins = int(60*dm)
1459
1460 ds = 60*dm - mins
1461 secs = int(60*ds)
1462

◆ diffusefraction()

subroutine beers_module::diffusefraction ( real(kind(1d0)), intent(in) radg,
real(kind(1d0)), intent(in) altitude,
real(kind(1d0)), intent(in) kt,
real(kind(1d0)), intent(in) ta,
real(kind(1d0)), intent(in) rh,
real(kind(1d0)), intent(out) radi,
real(kind(1d0)), intent(out) radd )

Definition at line 777 of file suews_phys_beers.f95.

778 IMPLICIT NONE
779
780 REAL(KIND(1D0)), INTENT(in) :: radG
781 REAL(KIND(1D0)), INTENT(in) :: altitude
782 REAL(KIND(1D0)), INTENT(in) :: Kt
783 REAL(KIND(1D0)), INTENT(in) :: Ta
784 REAL(KIND(1D0)), INTENT(in) :: RH
785 REAL(KIND(1D0)), INTENT(out) :: radD ! direct radiation
786 REAL(KIND(1D0)), INTENT(out) :: radI ! diffusive radiation
787 REAL(KIND(1D0)) :: alfa
788
789 alfa = altitude*deg2rad
790
791 IF (ta <= -99 .OR. rh <= -99) THEN !.or. isnan(Ta) .or. isnan(RH)) then
792 IF (kt <= 0.3) THEN
793 radd = radg*(1.020 - 0.248*kt)
794 ELSE IF (kt > 0.3 .AND. kt < 0.78) THEN
795 radd = radg*(1.45 - 1.67*kt)
796 ELSE IF (kt >= 0.78) THEN
797 radd = radg*0.147
798 END IF
799 ELSE
800 IF (kt <= 0.3) THEN
801 radd = radg*(1 - 0.232*kt + 0.0239*sin(alfa) - 0.000682*ta + 0.0195*(rh/100))
802 ELSE IF (kt > 0.3 .AND. kt < 0.78) THEN
803 radd = radg*(1.329 - 1.716*kt + 0.267*sin(alfa) - 0.00357*ta + 0.106*(rh/100))
804 ELSE IF (kt >= 0.78) THEN
805 radd = radg*(0.426*kt - 0.256*sin(alfa) + 0.00349*ta + 0.0734*(rh/100))
806 END IF
807 END IF
808 ! correction of radD
809 radd = max(0.d0, radd)
810 radd = min(radg, radd)
811
812 ! calculation of direct beam radiation
813 radi = (radg - radd)/(sin(alfa))
814
815 !! Corrections for low sun altitudes (20130307)
816 IF (altitude < 1 .AND. radi > radg) THEN
817 radi = radg
818 END IF
819

References allocatearray::deg2rad.

Referenced by beers_cal_main(), and tsurfbeers().

Here is the caller graph for this function:

◆ hwtosvf_ground()

real(kind(1d0)) function beers_module::hwtosvf_ground ( real(kind(1d0)), intent(in) hw)

Definition at line 457 of file suews_phys_beers.f95.

458 IMPLICIT NONE
459 REAL(KIND(1D0)), PARAMETER :: a = 0.5598
460 REAL(KIND(1D0)), PARAMETER :: b = -0.2485
461 REAL(KIND(1D0)), PARAMETER :: c = 0.4112
462 REAL(KIND(1D0)), PARAMETER :: d = -0.02528
463
464 REAL(KIND(1D0)), INTENT(in) :: hw
465 REAL(KIND(1D0)) :: svfGround
466
467 ! SvfGround: Parameterisation based on NYC data (500x500 meter grid)
468 svfground = a*exp(b*hw) + c*exp(d*hw)
469

Referenced by beers_cal_main().

Here is the caller graph for this function:

◆ hwtosvf_roof()

real(kind(1d0)) function beers_module::hwtosvf_roof ( real(kind(1d0)), intent(in) hw)

Definition at line 472 of file suews_phys_beers.f95.

473 IMPLICIT NONE
474 REAL(KIND(1D0)), PARAMETER :: e = 0.5572
475 REAL(KIND(1D0)), PARAMETER :: f = 0.0589
476 REAL(KIND(1D0)), PARAMETER :: g = 0.4143
477
478 REAL(KIND(1D0)), INTENT(in) :: hw
479 REAL(KIND(1D0)) :: svfRoof
480
481 ! SvfGround: Parameterisation based on NYC data (500x500 meter grid)
482 svfroof = e*exp(-f*hw) + g
483

Referenced by beers_cal_main().

Here is the caller graph for this function:

◆ issign()

subroutine beers_module::issign ( real(kind(1d0)) ix,
real(kind(1d0)) maxpos,
real(kind(1d0)) isignm )

Definition at line 1236 of file suews_phys_beers.f95.

1237 REAL(KIND(1D0)) IX, MAXPOS, ISIGNM
1238 isignm = 1.0
1239 IF (ix < 0 .OR. ix > maxpos) isignm = -1
1240 IF (ix == 0) isignm = 0
1241 RETURN

◆ kroof()

subroutine beers_module::kroof ( real(kind(1d0)), intent(in) radi,
real(kind(1d0)), intent(in) radd,
real(kind(1d0)), intent(in) radg,
real(kind(1d0)), intent(in) f_sh,
real(kind(1d0)), intent(in) altitude,
real(kind(1d0)), intent(in) svfr,
real(kind(1d0)), intent(in) svfveg,
real(kind(1d0)), intent(in) shadow,
real(kind(1d0)), intent(in) psi,
real(kind(1d0)), intent(in) alb_bldg,
real(kind(1d0)), intent(out) kdown )

Definition at line 378 of file suews_phys_beers.f95.

381
382 IMPLICIT NONE
383
384 REAL(KIND(1D0)), INTENT(in) :: radI
385 REAL(KIND(1D0)), INTENT(in) :: radG
386 REAL(KIND(1D0)), INTENT(in) :: radD
387 ! REAL(KIND(1D0)), intent(in) :: zen
388 REAL(KIND(1D0)), INTENT(in) :: altitude
389 REAL(KIND(1D0)), INTENT(in) :: svfr
390 REAL(KIND(1D0)), INTENT(in) :: psi
391 REAL(KIND(1D0)), INTENT(in) :: alb_bldg
392 REAL(KIND(1D0)), INTENT(in) :: shadow, svfveg
393 REAL(KIND(1D0)), INTENT(in) :: F_sh
394 REAL(KIND(1D0)), INTENT(out) :: Kdown
395
396 ! REAL(KIND(1D0)) :: svfalfa, tmp
397 REAL(KIND(1D0)) :: svfrbuveg
398
399 !Building height angle from svfr
400 ! tmp = 1 - (svfr + svfveg - 1)
401 ! if (tmp <= 0) tmp = 0.000000001 ! avoiding log(0)
402 ! svfalfa = ASIN(EXP(LOG(tmp)/2))
403 ! call cal_svfalfa(svfr, svfveg, svfalfa, tmp)
404
405 ! ! Fraction shadow on building walls based on sun altitude and svf
406 ! IF (altitude > 0) THEN
407 ! call cylindric_wedge(zen, svfalfa, F_sh)
408 ! if (F_sh < 0) F_sh = 0.0
409 ! if (F_sh > 0.5) F_sh = 0.5 !TODO #8 why there is no such restriction in python code?
410 ! else
411 ! F_sh = 1.
412 ! end if
413
414 svfrbuveg = svfr - (1.0 - svfveg)*(1.0 - psi)
415
416 kdown = radi*shadow*sin(altitude*deg2rad) &
417 + radd*svfrbuveg &
418 + alb_bldg*(1 - svfrbuveg)*(radg*(1 - f_sh) + radd*f_sh)
419

References allocatearray::deg2rad.

Referenced by beers_cal_main().

Here is the caller graph for this function:

◆ kvikt_veg()

subroutine beers_module::kvikt_veg ( real(kind(1d0)), intent(in) svf,
real(kind(1d0)), intent(in) svfveg,
real(kind(1d0)), intent(in) vikttot,
real(kind(1d0)), intent(out) viktveg,
real(kind(1d0)), intent(out) viktwall )

Definition at line 922 of file suews_phys_beers.f95.

924
925 IMPLICIT NONE
926 REAL(KIND(1D0)), INTENT(in) :: vikttot
927 REAL(KIND(1D0)), INTENT(in) :: svf
928 REAL(KIND(1D0)), INTENT(in) :: svfveg
929 REAL(KIND(1D0)), INTENT(out) :: viktveg, viktwall
930 REAL(KIND(1D0)) :: svfvegbu
931
932 !! Least
933 ! viktwall = (vikttot &
934 ! - (63.227*svf**6 - 161.51*svf**5 &
935 ! + 156.91*svf**4 - 70.424*svf**3 &
936 ! + 16.773*svf**2 - 0.4863*svf))/vikttot
937
938 viktwall = cal_vikt(svf, vikttot)
939
940 svfvegbu = (svfveg + svf - 1) ! Vegetation plus buildings
941 ! viktveg = (vikttot &
942 ! - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
943 ! + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
944 ! + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
945 viktveg = cal_vikt(svfvegbu, vikttot)
946 viktveg = viktveg - viktwall

References cal_vikt().

Referenced by kwalls().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ kwalls()

subroutine beers_module::kwalls ( real(kind(1d0)), intent(in) svf,
real(kind(1d0)), intent(in) svfveg,
real(kind(1d0)), intent(in) shadow,
real(kind(1d0)), intent(in) f_sh,
real(kind(1d0)), intent(in) radi,
real(kind(1d0)), intent(in) radg,
real(kind(1d0)), intent(in) radd,
real(kind(1d0)), intent(in) azimuth,
real(kind(1d0)), intent(in) altitude,
real(kind(1d0)), intent(in) psi,
real(kind(1d0)), intent(in) t,
real(kind(1d0)), intent(in) alb_ground,
real(kind(1d0)), intent(in) alb_bldg,
real(kind(1d0)), intent(out) keast,
real(kind(1d0)), intent(out) knorth,
real(kind(1d0)), intent(out) ksouth,
real(kind(1d0)), intent(out) kwest )

Definition at line 822 of file suews_phys_beers.f95.

826
827 IMPLICIT NONE
828
829 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
830 REAL(KIND(1D0)), INTENT(in) :: radI
831 REAL(KIND(1D0)), INTENT(in) :: radG
832 REAL(KIND(1D0)), INTENT(in) :: radD
833 REAL(KIND(1D0)), INTENT(in) :: azimuth
834 REAL(KIND(1D0)), INTENT(in) :: altitude
835 REAL(KIND(1D0)), INTENT(in) :: psi
836 REAL(KIND(1D0)), INTENT(in) :: t
837 REAL(KIND(1D0)), INTENT(in) :: alb_bldg
838 REAL(KIND(1D0)), INTENT(in) :: alb_ground
839
840 ! REAL(KIND(1d0)), intent(in) :: shadow, F_sh, svfalfa, svf, svfveg, svfaveg
841 REAL(KIND(1D0)), INTENT(in) :: shadow, F_sh, svf, svfveg
842 REAL(KIND(1D0)), INTENT(out) :: Keast, Knorth, Ksouth, Kwest
843
844 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW
845 REAL(KIND(1D0)) :: viktveg, viktwall
846 REAL(KIND(1D0)) :: KeastI, KsouthI, KwestI, KnorthI, Kuptowall
847 REAL(KIND(1D0)) :: KeastDG, KsouthDG, KwestDG, KnorthDG
848 REAL(KIND(1D0)) :: svfE, svfS, svfW, svfN, svfEveg, svfSveg, svfWveg, svfNveg, svfbuveg
849 REAL(KIND(1D0)) :: gvfalb, gvfalbnosh
850
851 ! Internal grids
852 REAL(KIND(1D0)) :: svfviktbuveg
853
854 ! ALLOCATE (svfviktbuveg(1, 1))
855 vikttot = 4.4897
856 azie = azimuth + t
857 azis = azimuth - 90 + t
858 aziw = azimuth - 180 + t
859 azin = azimuth - 270 + t
860
861 svfe = svf
862 svfs = svf
863 svfw = svf
864 svfn = svf
865 svfeveg = svfveg
866 svfsveg = svfveg
867 svfwveg = svfveg
868 svfnveg = svfveg
869
870 !!! Direct irradiance !!!
871 IF (azimuth > (360 - t) .OR. azimuth <= (180 - t)) THEN
872 keasti = radi*shadow*cos(altitude*deg2rad)*sin(azie*deg2rad)
873 ELSE
874 keasti = 0
875 END IF
876 IF (azimuth > (90 - t) .AND. azimuth <= (270 - t)) THEN
877 ksouthi = radi*shadow*cos(altitude*deg2rad)*sin(azis*deg2rad)
878 ELSE
879 ksouthi = 0
880 END IF
881 IF (azimuth > (180 - t) .AND. azimuth <= (360 - t)) THEN
882 kwesti = radi*shadow*cos(altitude*deg2rad)*sin(aziw*deg2rad)
883 ELSE
884 kwesti = 0
885 END IF
886 IF (azimuth <= (90 - t) .OR. azimuth > (270 - t)) THEN
887 knorthi = radi*shadow*cos(altitude*deg2rad)*sin(azin*deg2rad)
888 ELSE
889 knorthi = 0
890 END IF
891
892 !!! Diffuse and reflected radiation !!!
893 svfbuveg = svfs - (1.0 - svfsveg)*(1.0 - psi)
894 gvfalb = shadow*alb_ground ! albedo not defined TODO #3
895 gvfalbnosh = (1 - shadow)*alb_ground
896 kuptowall = (gvfalb*radi*sin(altitude*deg2rad)) &
897 + (radd*svfbuveg + alb_bldg*(1 - svfbuveg)*(radg*(1 - f_sh) + radd*f_sh))*gvfalbnosh
898
899 CALL kvikt_veg(svfe, svfeveg, vikttot, viktveg, viktwall)
900 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
901 keastdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
902 keast = keasti + keastdg
903
904 CALL kvikt_veg(svfs, svfsveg, vikttot, viktveg, viktwall)
905 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
906 ksouthdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
907 ksouth = ksouthi + ksouthdg
908
909 CALL kvikt_veg(svfw, svfwveg, vikttot, viktveg, viktwall)
910 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
911 kwestdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
912 kwest = kwesti + kwestdg
913
914 CALL kvikt_veg(svfn, svfnveg, vikttot, viktveg, viktwall)
915 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
916 knorthdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
917 knorth = knorthi + knorthdg
918
919 ! DEALLOCATE (svfviktbuveg)

References allocatearray::deg2rad, and kvikt_veg().

Referenced by beers_cal_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ leapyearcalc()

subroutine beers_module::leapyearcalc ( integer year_int,
integer nrodays )

Definition at line 1381 of file suews_phys_beers.f95.

1382
1383 IMPLICIT NONE
1384
1385 INTEGER :: nroDays, year_int
1386
1387 IF (mod(year_int, 100) /= 0 .AND. mod(year_int, 4) == 0) THEN
1388 nrodays = 366
1389 ELSEIF (mod(year_int, 400) == 0) THEN
1390 nrodays = 366
1391 ELSE
1392 nrodays = 365
1393 END IF

◆ lvikt_veg()

subroutine beers_module::lvikt_veg ( real(kind(1d0)), intent(in) isvf,
real(kind(1d0)), intent(in) isvfveg,
real(kind(1d0)), intent(in) isvfaveg,
real(kind(1d0)), intent(in) vikttot,
real(kind(1d0)), intent(out) viktveg,
real(kind(1d0)), intent(out) viktsky,
real(kind(1d0)), intent(out) viktrefl,
real(kind(1d0)), intent(out) viktwall )

Definition at line 1189 of file suews_phys_beers.f95.

1191
1192 IMPLICIT NONE
1193
1194 REAL(KIND(1D0)), INTENT(in) :: vikttot
1195 REAL(KIND(1D0)), INTENT(in) :: isvf
1196 REAL(KIND(1D0)), INTENT(in) :: isvfveg
1197 REAL(KIND(1D0)), INTENT(in) :: isvfaveg
1198 REAL(KIND(1D0)), INTENT(out) :: viktveg
1199 REAL(KIND(1D0)), INTENT(out) :: viktsky
1200 REAL(KIND(1D0)), INTENT(out) :: viktrefl
1201 REAL(KIND(1D0)), INTENT(out) :: viktwall
1202
1203 REAL(KIND(1D0)) :: viktonlywall
1204 REAL(KIND(1D0)) :: viktaveg
1205 REAL(KIND(1D0)) :: svfvegbu
1206
1207 viktonlywall = (vikttot - &
1208 (63.227*isvf**6 - 161.51*isvf**5 + 156.91*isvf**4 &
1209 - 70.424*isvf**3 + 16.773*isvf**2 - 0.4863*isvf))/vikttot
1210
1211 viktaveg = (vikttot &
1212 - (63.227*isvfaveg**6 - 161.51*isvfaveg**5 &
1213 + 156.91*isvfaveg**4 - 70.424*isvfaveg**3 &
1214 + 16.773*isvfaveg**2 - 0.4863*isvfaveg))/vikttot
1215
1216 viktwall = viktonlywall - viktaveg
1217
1218 svfvegbu = (isvfveg + isvf - 1) ! Vegetation plus buildings
1219 viktsky = (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1220 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1221 + 16.773*svfvegbu**2 - 0.4863*svfvegbu)/vikttot
1222 viktrefl = (vikttot &
1223 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1224 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1225 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1226 viktveg = (vikttot &
1227 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1228 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1229 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1230 viktveg = viktveg - viktwall
1231

Referenced by lwalls().

Here is the caller graph for this function:

◆ lwalls()

subroutine beers_module::lwalls ( real(kind(1d0)), intent(in) svf,
real(kind(1d0)), intent(in) svfveg,
real(kind(1d0)), intent(in) svfaveg,
real(kind(1d0)), intent(in) ldown2d,
real(kind(1d0)), intent(in) lup2d,
real(kind(1d0)), intent(in) altitude,
real(kind(1d0)), intent(in) ta,
real(kind(1d0)), intent(in) tw,
real(kind(1d0)), intent(in) sbc,
real(kind(1d0)), intent(in) emis_wall,
real(kind(1d0)), intent(in) emis_sky,
real(kind(1d0)), intent(in) t,
real(kind(1d0)), intent(in) ci,
real(kind(1d0)), intent(in) azimuth,
real(kind(1d0)), intent(in) ldown,
real(kind(1d0)), intent(in) svfalfa,
real(kind(1d0)), intent(in) f_sh_in,
real(kind(1d0)), intent(out) least,
real(kind(1d0)), intent(out) lnorth,
real(kind(1d0)), intent(out) lsouth,
real(kind(1d0)), intent(out) lwest )

Definition at line 961 of file suews_phys_beers.f95.

965 IMPLICIT NONE
966
967 REAL(KIND(1D0)), INTENT(in) :: altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, ldown
968
969 REAL(KIND(1D0)), INTENT(in) :: svfalfa, svf, svfveg, svfaveg, F_sh_in
970 REAL(KIND(1D0)), INTENT(in) :: Ldown2d, Lup2d
971 REAL(KIND(1D0)), INTENT(out) :: Least, Lnorth, Lsouth, Lwest
972
973 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW, c
974
975 REAL(KIND(1D0)) :: svfalfaE, svfalfaS, svfalfaW, svfalfaN
976 REAL(KIND(1D0)) :: alfaB, betaB, betasun
977 REAL(KIND(1D0)) :: Lground, Lrefl, Lsky, Lsky_allsky, Lveg, Lwallsh, Lwallsun
978 ! REAL(KIND(1d0)) :: viktonlywall, viktaveg, svfvegbu
979 ! REAL(KIND(1d0)) :: oneminussvfE, oneminussvfS, oneminussvfW, oneminussvfN
980 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
981 ! INTEGER, PARAMETER:: SOLWEIG_ldown = 0 !! force to 0, TS 13 Dec 2019
982
983 ! set as 1 for testing
984 INTEGER, PARAMETER :: SOLWEIG_ldown = 1 !! force to 0, TS 13 Dec 2019
985
986 REAL(KIND(1D0)) :: viktveg, viktsky, viktrefl, viktwall
987 REAL(KIND(1D0)) :: svfE, svfS, svfW, svfN, svfEveg, svfSveg, svfWveg, svfNveg
988 REAL(KIND(1D0)) :: svfEaveg, svfSaveg, svfWaveg, svfNaveg, F_sh
989
990 ! ALLOCATE (oneminussvfE(1, 1))
991 ! ALLOCATE (oneminussvfS(1, 1))
992 ! ALLOCATE (oneminussvfW(1, 1))
993 ! ALLOCATE (oneminussvfN(1, 1))
994 ! ALLOCATE (svfalfaE(1, 1))
995 ! ALLOCATE (svfalfaS(1, 1))
996 ! ALLOCATE (svfalfaW(1, 1))
997 ! ALLOCATE (svfalfaN(1, 1))
998 ! ALLOCATE (alfaB(1, 1))
999 ! ALLOCATE (betaB(1, 1))
1000 ! ALLOCATE (betasun(1, 1))
1001 ! ALLOCATE (Lground(1, 1))
1002 ! ALLOCATE (Lrefl(1, 1))
1003 ! ALLOCATE (Lsky(1, 1))
1004 ! ALLOCATE (Lsky_allsky(1, 1))
1005 ! ALLOCATE (Lveg(1, 1))
1006 ! ALLOCATE (Lwallsh(1, 1))
1007 ! ALLOCATE (Lwallsun(1, 1))
1008 ! ALLOCATE (viktonlywall(1, 1))
1009 ! ALLOCATE (viktaveg(1, 1))
1010 ! ALLOCATE (svfvegbu(1, 1))
1011
1012 svfe = svf
1013 svfs = svf
1014 svfw = svf
1015 svfn = svf
1016 svfeveg = svfveg
1017 svfsveg = svfveg
1018 svfwveg = svfveg
1019 svfnveg = svfveg
1020 svfeaveg = svfaveg
1021 svfsaveg = svfaveg
1022 svfwaveg = svfaveg
1023 svfnaveg = svfaveg
1024
1025 !Building height angle from svf
1026 ! oneminussvfE = 1.-svfE; WHERE (oneminussvfE <= 0) oneminussvfE = 0.000000001 ! avoiding log(0)
1027 ! oneminussvfS = 1.-svfS; WHERE (oneminussvfS <= 0) oneminussvfS = 0.000000001 ! avoiding log(0)
1028 ! oneminussvfW = 1.-svfW; WHERE (oneminussvfW <= 0) oneminussvfW = 0.000000001 ! avoiding log(0)
1029 ! oneminussvfN = 1.-svfN; WHERE (oneminussvfN <= 0) oneminussvfN = 0.000000001 ! avoiding log(0)
1030 ! svfalfaE = ASIN(EXP((LOG(oneminussvfE))/2))
1031 ! svfalfaS = ASIN(EXP((LOG(oneminussvfS))/2))
1032 ! svfalfaW = ASIN(EXP((LOG(oneminussvfW))/2))
1033 ! svfalfaN = ASIN(EXP((LOG(oneminussvfN))/2))
1034 svfalfae = svfalfa
1035 svfalfas = svfalfa
1036 svfalfaw = svfalfa
1037 svfalfan = svfalfa
1038
1039 vikttot = 4.4897
1040 aziw = azimuth + t
1041 azin = azimuth - 90 + t
1042 azie = azimuth - 180 + t
1043 azis = azimuth - 270 + t
1044
1045 f_sh = 2.*f_sh_in - 1. !(cylindric_wedge scaled 0-1 as only half hemisphere is seen)
1046
1047 IF (solweig_ldown == 1) THEN
1048 c = 1 - ci
1049 lsky_allsky = emis_sky*sbc*((ta + 273.15)**4)*(1 - c) + c*sbc*((ta + 273.15)**4)
1050 ELSE
1051 lsky_allsky = ldown
1052 END IF
1053
1054 !! Least
1055 CALL lvikt_veg(svfe, svfeveg, svfeaveg, vikttot, &
1056 viktveg, viktsky, viktrefl, viktwall)
1057
1058 IF (altitude > 0) THEN ! daytime
1059 alfab = atan(svfalfae)
1060 betab = atan(tan((svfalfae)*f_sh))
1061 betasun = ((alfab - betab)/2) + betab
1062 IF (azimuth > (180 - t) .AND. azimuth <= (360 - t)) THEN
1063 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azie*(pi/180)))**4)* &
1064 viktwall*(1 - f_sh)*cos(betasun)*0.5
1065 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1066 ELSE
1067 lwallsun = 0
1068 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1069 END IF
1070 ELSE !nighttime
1071 lwallsun = 0
1072 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1073 END IF
1074 lsky = ((svfe + svfeveg - 1)*lsky_allsky)*viktsky*0.5
1075 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1076 lground = lup2d*0.5
1077 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1078 least = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1079
1080 !! Lsouth
1081 CALL lvikt_veg(svfs, svfsveg, svfsaveg, vikttot, &
1082 viktveg, viktsky, viktrefl, viktwall)
1083
1084 IF (altitude > 0) THEN ! daytime
1085 alfab = atan(svfalfas)
1086 betab = atan(tan((svfalfas)*f_sh))
1087 betasun = ((alfab - betab)/2) + betab
1088 IF (azimuth <= (90 - t) .OR. azimuth > (270 - t)) THEN
1089 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azis*(pi/180)))**4)* &
1090 viktwall*(1 - f_sh)*cos(betasun)*0.5
1091 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1092 ELSE
1093 lwallsun = 0
1094 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1095 END IF
1096 ELSE !nighttime
1097 lwallsun = 0
1098 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1099 END IF
1100 lsky = ((svfs + svfsveg - 1)*lsky_allsky)*viktsky*0.5
1101 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1102 lground = lup2d*0.5
1103 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1104 lsouth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1105
1106 !! Lwest
1107 CALL lvikt_veg(svfw, svfwveg, svfwaveg, vikttot, &
1108 viktveg, viktsky, viktrefl, viktwall)
1109
1110 IF (altitude > 0) THEN ! daytime
1111 alfab = atan(svfalfaw)
1112 betab = atan(tan((svfalfaw)*f_sh))
1113 betasun = ((alfab - betab)/2) + betab
1114 IF (azimuth > (360 - t) .OR. azimuth <= (180 - t)) THEN
1115 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(aziw*(pi/180)))**4)* &
1116 viktwall*(1 - f_sh)*cos(betasun)*0.5
1117 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1118 ELSE
1119 lwallsun = 0
1120 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1121 END IF
1122 ELSE !nighttime
1123 lwallsun = 0
1124 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1125 END IF
1126 lsky = ((svfw + svfwveg - 1)*lsky_allsky)*viktsky*0.5
1127 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1128 lground = lup2d*0.5
1129 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1130 lwest = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1131
1132 !! Lnorth
1133 CALL lvikt_veg(svfn, svfnveg, svfnaveg, vikttot, &
1134 viktveg, viktsky, viktrefl, viktwall)
1135
1136 IF (altitude > 0) THEN ! daytime
1137 alfab = atan(svfalfan)
1138 betab = atan(tan((svfalfan)*f_sh))
1139 betasun = ((alfab - betab)/2) + betab
1140 IF (azimuth > (90 - t) .AND. azimuth <= (270 - t)) THEN
1141 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azin*(pi/180)))**4)* &
1142 viktwall*(1 - f_sh)*cos(betasun)*0.5
1143 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1144 ELSE
1145 lwallsun = 0
1146 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1147 END IF
1148 ELSE !nighttime
1149 lwallsun = 0
1150 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1151 END IF
1152 lsky = ((svfn + svfnveg - 1)*lsky_allsky)*viktsky*0.5
1153 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1154 lground = lup2d*0.5
1155 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1156 lnorth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1157
1158 ! DEALLOCATE (svfE)
1159 ! DEALLOCATE (svfS)
1160 ! DEALLOCATE (svfW)
1161 ! DEALLOCATE (svfN)
1162 ! DEALLOCATE (svfEveg)
1163 ! DEALLOCATE (svfSveg)
1164 ! DEALLOCATE (svfWveg)
1165 ! DEALLOCATE (svfNveg)
1166 ! DEALLOCATE (svfEaveg)
1167 ! DEALLOCATE (svfSaveg)
1168 ! DEALLOCATE (svfWaveg)
1169 ! DEALLOCATE (svfNaveg)
1170 ! DEALLOCATE (svfalfaE)
1171 ! DEALLOCATE (svfalfaS)
1172 ! DEALLOCATE (svfalfaW)
1173 ! DEALLOCATE (svfalfaN)
1174 ! DEALLOCATE (alfaB)
1175 ! DEALLOCATE (betaB)
1176 ! DEALLOCATE (betasun)
1177 ! DEALLOCATE (Lground)
1178 ! DEALLOCATE (Lrefl)
1179 ! DEALLOCATE (Lsky)
1180 ! DEALLOCATE (Lsky_allsky)
1181 ! DEALLOCATE (Lveg)
1182 ! DEALLOCATE (Lwallsh)
1183 ! DEALLOCATE (Lwallsun)
1184 ! DEALLOCATE (viktonlywall)
1185 ! DEALLOCATE (viktaveg)
1186

References lvikt_veg(), and pi.

Referenced by beers_cal_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ month2day()

subroutine beers_module::month2day ( integer mon,
integer ne,
integer k,
integer b )

Definition at line 1348 of file suews_phys_beers.f95.

1349 IMPLICIT NONE
1350 INTEGER :: mon, ne, k, b
1351
1352 IF (mon == 1) THEN
1353 ne = 32 - b
1354 ELSE IF (mon == 2) THEN
1355 ne = 60 + k - b
1356 ELSE IF (mon == 3) THEN
1357 ne = 91 + k - b
1358 ELSE IF (mon == 4) THEN
1359 ne = 121 + k - b
1360 ELSE IF (mon == 5) THEN
1361 ne = 152 + k - b
1362 ELSE IF (mon == 6) THEN
1363 ne = 182 + k - b
1364 ELSE IF (mon == 7) THEN
1365 ne = 213 + k - b
1366 ELSE IF (mon == 8) THEN
1367 ne = 244 + k - b
1368 !**********PAGE 151 STARTS HERE**************
1369 ELSE IF (mon == 9) THEN
1370 ne = 274 + k - b
1371 ELSE IF (mon == 10) THEN
1372 ne = 305 + k - b
1373 ELSE IF (mon == 11) THEN
1374 ne = 335 + k - b
1375 ELSE IF (mon == 12) THEN
1376 ne = 366 + k - b
1377 END IF

◆ shadowgroundkusaka()

subroutine beers_module::shadowgroundkusaka ( real(kind(1d0)), intent(in) hw,
real(kind(1d0)), intent(in) azimuth,
real(kind(1d0)), intent(in) zen,
real(kind(1d0)), intent(out) shadowground,
real(kind(1d0)), intent(out) shadowwalls )

Definition at line 544 of file suews_phys_beers.f95.

545 ! Translated from Python! Should be checked by fortran person! FL
546 ! This function calculates sunlit fractions on ground and walls for an infinite long, rotating canyon (Kusaka)
547 ! (hw, azimuth, zen)
548 !DEG2RAD = np.pi / 180.
549
550 IMPLICIT NONE
551
552 REAL(KIND(1D0)), INTENT(in) :: HW, azimuth, zen
553 REAL(KIND(1D0)), INTENT(out) :: shadowground, shadowwalls
554 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
555 REAL(KIND(1D0)) :: ROOF_HEIGHT, ROAD_WIDTH, THETA_Z, THETA_S, THETA_S180, ndir, wx
556 REAL(KIND(1D0)) :: LshadowRoad(1:8), Wallsun(1:8)
557 INTEGER :: idir
558
559 roof_height = hw
560 road_width = 1.
561 theta_z = zen ! Solar Zenith Angle (radians)
562 theta_s = azimuth*deg2rad ! Solar azimuth Angle (radians)
563 theta_s180 = abs(pi - theta_s) ! deviation from south acc. Kusaka 2001
564 ndir = 8.
565 !LshadowRoad = np.zeros(int(ndir))
566 !Wallsun = np.zeros(int(ndir))
567
568 DO idir = 1, 8 ! Looping through 8 canyon directions. In Python 1 to 9. Here 1 to 8 Correct???
569 wx = 1./sin(abs((float(idir)*pi/ndir) - theta_s180)) ! sun azimuth length in canyon
570 lshadowroad(idir) = (roof_height/road_width)*tan(theta_z)*sin(abs((float(idir)*pi/ndir) - theta_s180))
571 wallsun(idir) = (wx/tan(theta_z))/(roof_height/road_width)
572 END DO
573
574 ! shadow fraction ground
575 WHERE (lshadowroad >= 1.) lshadowroad = 1.
576 lshadowroad = -1*lshadowroad + 1. ! Getting sunlit instead of shadow
577 shadowground = sum(lshadowroad)/ndir
578
579 ! shadow fraction walls
580 WHERE (wallsun >= 1.) wallsun = 1.
581 shadowwalls = (sum(wallsun)/ndir)/2. ! two walls which one is always in shadow
582

References allocatearray::deg2rad, and pi.

Referenced by beers_cal_main().

Here is the caller graph for this function:

◆ suews_cal_dectime()

subroutine beers_module::suews_cal_dectime ( integer, intent(in) id,
integer, intent(in) it,
integer, intent(in) imin,
integer, intent(in) isec,
real(kind(1d0)), intent(out) dectime )

Definition at line 1513 of file suews_phys_beers.f95.

1516 IMPLICIT NONE
1517 INTEGER, INTENT(in) :: id, it, imin, isec
1518
1519 REAL(KIND(1D0)), INTENT(out) :: dectime ! nsh in type real
1520
1521 dectime = real(id - 1, kind(1d0)) &
1522 + real(it, kind(1d0))/24 &
1523 + real(imin, kind(1d0))/(60*24) &
1524 + real(isec, kind(1d0))/(60*60*24)
1525

Referenced by suews_driver::suews_cal_main(), and suews_driver::suews_cal_main_dts().

Here is the caller graph for this function:

◆ suews_cal_dls()

subroutine beers_module::suews_cal_dls ( integer, intent(in) id,
integer, intent(in) startdls,
integer, intent(in) enddls,
integer, intent(out) dls )

Definition at line 1569 of file suews_phys_beers.f95.

1572 IMPLICIT NONE
1573
1574 INTEGER, INTENT(in) :: id, startDLS, endDLS
1575 INTEGER, INTENT(out) :: DLS
1576
1577 dls = 0
1578 IF (id > startdls .AND. id < enddls) dls = 1
1579

Referenced by suews_driver::suews_cal_main(), and suews_driver::suews_cal_main_dts().

Here is the caller graph for this function:

◆ suews_cal_tstep()

subroutine beers_module::suews_cal_tstep ( integer, intent(in) tstep,
integer, intent(out) nsh,
real(kind(1d0)), intent(out) nsh_real,
real(kind(1d0)), intent(out) tstep_real )

Definition at line 1529 of file suews_phys_beers.f95.

1532 IMPLICIT NONE
1533 INTEGER, INTENT(in) :: tstep ! number of timesteps per hour
1534 ! values that are derived from tstep
1535 INTEGER, INTENT(out) :: nsh ! number of timesteps per hour
1536 REAL(KIND(1D0)), INTENT(out) :: nsh_real ! nsh in type real
1537 REAL(KIND(1D0)), INTENT(out) :: tstep_real ! tstep in type real
1538 nsh = 3600/tstep
1539 nsh_real = nsh*1.0
1540 tstep_real = tstep*1.0
1541

Referenced by suews_driver::suews_cal_main(), and suews_driver::suews_cal_main_dts().

Here is the caller graph for this function:

◆ suews_cal_weekday()

subroutine beers_module::suews_cal_weekday ( integer, intent(in) iy,
integer, intent(in) id,
real(kind(1d0)), intent(in) lat,
integer, dimension(3), intent(out) dayofweek_id )

Definition at line 1544 of file suews_phys_beers.f95.

1547 IMPLICIT NONE
1548
1549 INTEGER, INTENT(in) :: iy ! year
1550 INTEGER, INTENT(in) :: id ! day of year
1551 REAL(KIND(1D0)), INTENT(in) :: lat
1552
1553 INTEGER, DIMENSION(3), INTENT(OUT) :: dayofWeek_id
1554
1555 INTEGER :: wd
1556 INTEGER :: mb
1557 INTEGER :: date
1558 INTEGER :: seas
1559
1560 CALL day2month(id, mb, date, seas, iy, lat) !Calculate real date from doy
1561 CALL day_of_week(date, mb, iy, wd) !Calculate weekday (1=Sun, ..., 7=Sat)
1562
1563 dayofweek_id(1) = wd !Day of week
1564 dayofweek_id(2) = mb !Month
1565 dayofweek_id(3) = seas !Season
1566
day_of_week
subroutine day_of_week(date, month, year, dow)
Definition suews_util_time.f95:174
day2month
subroutine day2month(b, mb, md, seas, year, latitude)
Definition suews_util_time.f95:24

References day2month(), and day_of_week().

Referenced by suews_driver::suews_cal_main(), and suews_driver::suews_cal_main_dts().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sun_distance()

subroutine beers_module::sun_distance ( integer jday,
real(kind(1d0)) d )

Definition at line 692 of file suews_phys_beers.f95.

693
694 ! Calculates solar irradiance variation based on mean earth sun distance
695 ! with day of year as input.
696 ! Partridge and Platt, 1975
697
698 INTEGER :: jday
699 REAL(KIND(1D0)) :: b, D
700
701 b = 2*3.141592654*jday/365
702 d = sqrt(1.00011 + 0.034221*cos(b) + 0.001280*sin(b) + 0.000719*cos(2*b) + 0.000077*sin(2*b))
703

Referenced by clearnessindex_2013b().

Here is the caller graph for this function:

◆ tsurfbeers()

subroutine beers_module::tsurfbeers ( integer, intent(in) iy,
real(kind(1d0)), intent(in) ta,
real(kind(1d0)), intent(in) rh,
real(kind(1d0)), intent(in) radi,
real(kind(1d0)), intent(in) i0,
real(kind(1d0)), intent(in) dectime,
real(kind(1d0)), intent(in) snup,
real(kind(1d0)), intent(in) altitude,
real(kind(1d0)), intent(in) zen,
real(kind(1d0)), intent(in) timezone,
real(kind(1d0)), intent(in) lat,
real(kind(1d0)), intent(in) lng,
real(kind(1d0)), intent(in) alt,
real(kind(1d0)), intent(out) tg,
real(kind(1d0)), intent(out) tgwall,
real(kind(1d0)), intent(out) altmax )

Definition at line 486 of file suews_phys_beers.f95.

488 ! Ground surface temperature wave from SOLWEIG
489
490 IMPLICIT NONE
491
492 INTEGER, INTENT(in) :: iy
493 REAL(KIND(1D0)), INTENT(in) :: Ta, RH, radI, I0
494 REAL(KIND(1D0)), INTENT(in) :: dectime, SNUP, altitude, zen
495 REAL(KIND(1D0)), INTENT(in) :: timezone, lat, lng, alt
496 REAL(KIND(1D0)), INTENT(out) :: Tg, Tgwall, altmax
497 REAL(KIND(1D0)), PARAMETER :: TgK = 0.37
498 REAL(KIND(1D0)), PARAMETER :: Tstart = 3.41
499 REAL(KIND(1D0)), PARAMETER :: TgK_wall = 0.37
500 REAL(KIND(1D0)), PARAMETER :: Tstart_wall = -3.41
501 REAL(KIND(1D0)), PARAMETER :: TmaxLST = 15.
502 REAL(KIND(1D0)), PARAMETER :: TmaxLST_wall = 15.
503
504 REAL(KIND(1D0)) :: Ktc, notU, Tgamp, Tgampwall, radI0, corr, CI_Tg
505 REAL(KIND(1D0)) :: fifteen, sunmaximum, zen_sunmax, dectimemax, azimuth
506
507 ! finding daily maximum solar altitude
508 fifteen = 0.
509 sunmaximum = -90.
510 zen_sunmax = 90.
511 DO WHILE (sunmaximum <= 90.-zen_sunmax)
512 sunmaximum = 90.-zen_sunmax
513 fifteen = fifteen + 15./1440.
514 dectimemax = floor(dectime) + 10/24.+fifteen
515 CALL narp_cal_sunposition( &
516 REAL(iy, KIND(1D0)), & !input:
517 dectimemax, & ! sun position at middle of timestep before
518 timezone, lat, lng, alt, &
519 azimuth, zen_sunmax) !output:
520 END DO
521
522 altmax = sunmaximum
523
524 tgamp = (tgk*altmax - tstart) + tstart
525 tgampwall = (tgk_wall*altmax - (tstart_wall)) + (tstart_wall)
526 tg = tgamp*sin((((dectime - floor(dectime)) - snup/24)/(tmaxlst/24 - snup/24))*pi/2) + tstart
527 tgwall = tgampwall*sin((((dectime - floor(dectime)) - snup/24)/(tmaxlst_wall/24 - snup/24))*pi/2) + (tstart_wall)
528
529 ktc = 1.0
530 CALL diffusefraction(i0, altitude, ktc, ta, rh, radi0, notu)
531 corr = 0.1473*log(90.-(zen*rad2deg)) + 0.3454 ! 20070329 temporary correction of latitude from Lindberg et al. 2008
532 ci_tg = (radi/radi0) + (1.-corr)
533
534 !TODO: FIX THIS in fortran syntax?? if CI_Tg<inf then CI_Tg=1
535 IF (ci_tg > 1) ci_tg = 1
536 tg = tg*ci_tg
537 tgwall = tgwall*ci_tg
538
539 IF (tgwall < 0) tgwall = 0
540 IF (tg < 0) tg = 0
541

References diffusefraction(), narp_module::narp_cal_sunposition(), pi, and allocatearray::rad2deg.

Referenced by beers_cal_main().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ deg2rad

real(kind(1d0)), parameter beers_module::deg2rad = pi/180

Definition at line 32 of file suews_phys_beers.f95.

32 REAL(KIND(1D0)), PARAMETER :: DEG2RAD = pi/180

◆ pi

real(kind(1d0)), parameter beers_module::pi = ATAN(1.)*4

Definition at line 31 of file suews_phys_beers.f95.

31 REAL(KIND(1D0)), PARAMETER :: pi = atan(1.)*4

Referenced by cal_ci_latenight(), clearnessindex_2013b(), cylindric_wedge(), daylen(), lwalls(), shadowgroundkusaka(), and tsurfbeers().

◆ rad2deg

real(kind(1d0)), parameter beers_module::rad2deg = 1/DEG2RAD

Definition at line 33 of file suews_phys_beers.f95.

33 REAL(KIND(1D0)), PARAMETER :: RAD2deg = 1/deg2rad
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0