SUEWS API Site
Documentation of SUEWS source code
suews_phys_solweig.f95
Go to the documentation of this file.
1MODULE solweig_module
2 ! USE data_in
3 ! USE gis_data
4 ! USE time
5 USE allocatearray, ONLY: ncolumnsdataoutsolweig, deg2rad, rad2deg
6 ! USE InitialCond
7 ! USE sues_data
8 USE defaultnotused, ONLY: notused, notusedi
9
10 IMPLICIT NONE
11
12 ! REAL(KIND(1d0)) :: timestepdec !time step in decimal time
13 ! REAL(KIND(1d0)) :: CIlatenight
14 ! REAL(KIND(1d0)) :: timeadd
15 ! REAL(KIND(1d0)) :: firstdaytime ! if new day starts, =1 else =0
16 ! REAL(KIND(1d0)) :: Fside ! fraction of a person seen from each cardinal point
17 ! REAL(KIND(1d0)) :: Fup ! fraction of a person seen from down and up
18 ! REAL(KIND(1d0)) :: scale
19 ! REAL(KIND(1d0)) :: amaxvalue
20 ! REAL(KIND(1d0)) :: trans
21 ! REAL(KIND(1d0)) :: transperLAI
22 ! REAL(KIND(1d0)) :: xllcorner
23 ! REAL(KIND(1d0)) :: yllcorner
24 ! REAL(KIND(1d0)) :: NoData
25 ! REAL(KIND(1d0)) :: cellsize
26 ! INTEGER :: SolweigCount
27
28 ! INTEGER, PARAMETER:: SOLWEIG_ldown = 0 !! force not to use SOLWEIG based Ldown calculations, TS 13 Dec 2019
29
30 ! INTEGER, parameter :: 1 = 1 ! number of rows and cols of grid
31 ! INTEGER, parameter :: 1 = 1 ! number of rows and cols of grid
32
33 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: a, sh, vbshvegsh, vegsh
34 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: bush, vegdem, vegdem2, tempgrid
35 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: buildings, svf, svfE, svfS, svfW, svfN
36 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: svfveg, svfEveg, svfSveg, svfWveg, svfNveg
37 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg, last
38 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: Kdown2d, Keast, Knorth, Ksouth, Kup2d, Kwest
39 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: Ldown2d, Least, Lnorth, Lsouth, Lup2d, Lwest
40 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: gvf, Tmrt, shadow, Sstr, F_sh, sunwall
41 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: svfalfa, sos, Tgmap1
42 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: viktveg, viktsky, viktrefl, viktwall, savegrid
43CONTAINS
44 ! This is the core function of the SOLWEIG model
45 ! 2013-10-27
46 ! Fredrik Lindberg, fredrikl@gvc.gu.se
47 ! Gothenburg Urban Climate Group
48 ! Gothenburg University
49 ! Last modified:
50 ! LJ 27 Jan 2016 Renoval of tabs and fixing int-real conversions
51 ! HCW 02 Dec 2014 DEG2RAD and RAD2DEG commented out as now defined in AllocateArray
52
53 SUBROUTINE solweig_cal_main( &
54 id, &
55 it, &
56 dectime, &
57 lamdaP, &
58 lamdaF, &
59 avkdn, &
60 ldown, &
61 Temp_C, &
62 avrh, &
63 Press_hPa, &
64 Tg, &
65 lat, &
66 zenith_deg, &
67 azimuth, &
68 scale, &
69 alb_ground, &
70 alb_bldg, &
71 emis_ground, &
72 emis_wall, &
73 heightgravity, &
74 dataOutLineSOLWEIG) ! output
75
76 ! USE matsize
77 ! USE solweig_module
78
79 IMPLICIT NONE
80
81 INTEGER, INTENT(in) :: id
82 INTEGER, INTENT(in) :: it
83 REAL(KIND(1D0)), INTENT(in) :: lamdaP ! plan area fraction
84 REAL(KIND(1D0)), INTENT(in) :: lamdaF ! frontal area fraction
85
86 ! REAL(KIND(1d0)), intent(in) ::lai_id_dectr
87 ! REAL(KIND(1d0)), intent(in) ::LAImax_dectr
88 REAL(KIND(1D0)), INTENT(in) :: Press_hPa
89 REAL(KIND(1D0)), INTENT(in) :: Temp_C
90 REAL(KIND(1D0)), INTENT(in) :: avrh
91 REAL(KIND(1D0)), INTENT(in) :: avkdn
92 REAL(KIND(1D0)), INTENT(in) :: ldown
93 REAL(KIND(1D0)), INTENT(in) :: Tg
94 ! REAL(KIND(1d0)), intent(in) ::kdiff
95 ! REAL(KIND(1d0)), intent(in) ::kdir
96 REAL(KIND(1D0)), INTENT(in) :: lat
97 REAL(KIND(1D0)), INTENT(in) :: zenith_deg
98 REAL(KIND(1D0)), INTENT(in) :: azimuth
99 REAL(KIND(1D0)), INTENT(in) :: dectime
100 REAL(KIND(1D0)), INTENT(in) :: scale ! what is this?
101
102 REAL(KIND(1D0)), PARAMETER :: absL = 0.97 ! Absorption coefficient of longwave radiation of a person
103 REAL(KIND(1D0)), PARAMETER :: absK = 0.7 ! Absorption coefficient of shortwave radiation of a person
104 REAL(KIND(1D0)), INTENT(in) :: heightgravity ! Centre of gravity for a standing person
105 ! REAL(KIND(1D0)), intent(in)::TransMin ! Tranmissivity of K through decidious vegetation (leaf on)
106 REAL(KIND(1D0)), INTENT(in) :: alb_bldg ! Tranmissivity of K through decidious vegetation (leaf off)
107 REAL(KIND(1D0)), INTENT(in) :: alb_ground ! Tranmissivity of K through decidious vegetation (leaf off)
108 REAL(KIND(1D0)), INTENT(in) :: emis_wall ! Tranmissivity of K through decidious vegetation (leaf off)
109 REAL(KIND(1D0)), INTENT(in) :: emis_ground ! Tranmissivity of K through decidious vegetation (leaf off)
110
111 REAL(KIND(1D0)), DIMENSION(ncolumnsdataOutSOLWEIG - 5), INTENT(OUT) :: dataOutLineSOLWEIG
112
113 REAL(KIND(1D0)) :: t, Tstart, height, psi !,timezone,lat,lng,alt,amaxvalue
114 REAL(KIND(1D0)) :: altitude, zen !scale,azimuth,zenith
115 REAL(KIND(1D0)) :: CI, c, I0, Kt, Tw, weight1
116 REAL(KIND(1D0)) :: Ta, RH, P, radG, radD, radI !,idectime,tdectime!dectime,
117 REAL(KIND(1D0)) :: I0et, CIuncorr !,lati
118 REAL(KIND(1D0)) :: SNDN, SNUP, DEC, DAYL !,timestepdec,YEAR
119 REAL(KIND(1D0)) :: msteg, emis_sky, ea
120 ! REAL(KIND(1d0)),intent(in) ::lai_id
121 INTEGER :: DOY, hour, first, second, j !,ith!onlyglobal,usevegdem,x,y,i
122 REAL(KIND(1D0)) :: timeadd
123 REAL(KIND(1D0)) :: firstdaytime ! if new day starts, =1 else =0
124 ! REAL(KIND(1d0)) :: timestepdec !time step in decimal time
125 REAL(KIND(1D0)) :: CIlatenight
126 REAL(KIND(1D0)) :: Fside ! fraction of a person seen from each cardinal point
127 REAL(KIND(1D0)) :: Fup ! fraction of a person seen from down and up
128 REAL(KIND(1D0)) :: HW ! building height to width ratio
129
130 REAL(KIND(1D0)), DIMENSION(1, 1) :: svfalfa, sos, Tgmap1
131 REAL(KIND(1D0)), DIMENSION(1, 1) :: gvf !Ground View Factors (GVF)
132 REAL(KIND(1D0)), DIMENSION(1, 1) :: Tmrt, shadow, Sstr, F_sh
133 REAL(KIND(1D0)), DIMENSION(1, 1) :: vegsh
134 REAL(KIND(1D0)), DIMENSION(1, 1) :: buildings, svf
135 REAL(KIND(1D0)), DIMENSION(1, 1) :: svfveg
136 REAL(KIND(1D0)), DIMENSION(1, 1) :: svfaveg
137 REAL(KIND(1D0)), DIMENSION(1, 1) :: Kdown2d, Keast, Knorth, Ksouth, Kup2d, Kwest
138 REAL(KIND(1D0)), DIMENSION(1, 1) :: Ldown2d, Least, Lnorth, Lsouth, Lup2d, Lwest
139
140 ! Internal grids
141 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: tmp, Knight, svf_blgd_veg, Tgmap0 !,Tgmap
142 !Search directions for Ground View Factors (GVF)
143 REAL(KIND(1D0)), PARAMETER :: azimuthA(1:18) = [(j*(360.0/18.0), j=0, 17)]
144 ! temporary parameters and variables for testing
145 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
146 REAL(KIND(1D0)), PARAMETER :: SBC = 5.67051e-8
147 !REAL(KIND(1D0)),PARAMETER :: DEG2RAD=0.017453292,RAD2DEG=57.29577951 !Now defined in AllocateArray HCW 02 Dec 2014
148
149 INTEGER, PARAMETER :: onlyglobal = 1 !! force to calculate direct and diffuse components, TS 13 Dec 2019
150 INTEGER, PARAMETER :: usevegdem = 0 !! force not to use vegetation DEM based calculations, TS 13 Dec 2019
151 INTEGER, PARAMETER :: row = 1 !! force to 1, TS 13 Dec 2019
152 INTEGER, PARAMETER :: col = 1 !! force to 1, TS 13 Dec 2019
153 INTEGER, PARAMETER :: Posture = 1 !! force to 1, TS 13 Dec 2019
154 INTEGER, PARAMETER :: SOLWEIG_ldown = 0 !! force to 0, TS 13 Dec 2019
155
156 ! INTEGER:: firstTimeofDay = 0 !!Needs updating for new model timestep
157
158!!!!!! Begin program !!!!!!
159 ! internal grids
160 ALLOCATE (tmp(1, 1))
161 ALLOCATE (knight(1, 1))
162 ALLOCATE (tgmap0(1, 1))
163 ALLOCATE (svf_blgd_veg(1, 1))
164 !allocate(Tgmap0(1,1))
165
166 ! initialise this as ONE
167 cilatenight = 1
168
169 IF (posture == 1) THEN
170 fside = 0.22
171 fup = 0.06
172 ELSE
173 fside = 0.1666666667
174 fup = 0.166666667
175 END IF
176
177 ! These variables should change name in the future...
178 p = press_hpa
179 ta = temp_c
180 rh = avrh
181 radg = avkdn
182 doy = int(id)
183 hour = int(it)
184 height = heightgravity
185
186 ! TODO: need to check this
187 psi = 0
188
189 ! alb_bldg = alb(2) ! taken from Bldg (FunctionalTypes)
190 ! alb_ground = alb(1) ! taken from Paved (FunctionalTypes)
191 ! emis_wall = emis(2) ! taken from Bldg (FunctionalTypes)
192 ! emis_ground = emis(1) ! taken from Paved (FunctionalTypes)
193 ! radD = kdiff
194 ! radI = kdir
195
196 ! Transmissivity of shortwave radiation through vegetation based on decid LAI
197 ! trans = TransMin + (LAImax_dectr - LAI_id_dectr)*transperLAI
198 ! IF (it == firstTimeofDay) THEN
199 ! trans = TransMin + (LAImax_dectr - LAI_id_dectr)*transperLAI
200 ! ENDIF
201
202 ! Radiation sensor setup offset in degrees
203 t = 0
204
205 !Surface temperature difference at sunrise
206 tstart = 3.41
207
208 !Initialization of maps
209 knight = 0.0
210
211 hw = cal_ratio_height2width(lamdap, lamdaf)
212 ! parameterisation using NYC building data
213 svf = hwtosvf_ground(hw)
214
215 ! svfveg: SVF based on vegetation blocking the sky (1 = no vegetation)
216 svfveg = 1
217
218 ! TODO: what are these?
219 svfaveg = 1
220 buildings = 1
221
222 tmp = 1 - (svf + svfveg - 1)
223 WHERE (tmp <= 0) tmp = 0.000000001 ! avoiding log(0)
224 svfalfa = asin(exp(log(tmp)/2))
225
226 !Parameterization for Lup
227 first = int(anint(height)) !Radiative surface influence, Rule of thumb by Schmid et al. (1990).
228 IF (first == 0) THEN
229 first = 1
230 END IF
231 second = int(anint(height*20))
232
233 ! SVF combines for buildings and vegetation
234 svf_blgd_veg = (svf - (1 - svfveg)*(1 - psi))
235
236 ! Sun position related things
237 CALL daylen(doy, lat, dayl, dec, sndn, snup)
238 zen = zenith_deg*deg2rad
239 altitude = 90 - zenith_deg
240
241 !Determination of clear-sky emissivity from Prata (1996)
242 ea = 6.107*10**((7.5*ta)/(237.3 + ta))*(rh/100) !Vapor pressure
243 msteg = 46.5*(ea/(ta + 273.15))
244 emis_sky = (1 - (1 + msteg)*exp(-((1.2 + 3.0*msteg)**0.5))) - 0.04
245
246!!! DAYTIME !!!
247 IF (altitude > 0) THEN
248
249 !Clearness Index on Earth's surface after Crawford and Dunchon (1999) with a correction
250 !factor for low sun elevations after Lindberg et al. (2008)
251 CALL clearnessindex_2013b(zen, doy, ta, rh/100, radg, lat, p, i0, ci, kt, i0et, ciuncorr)
252 IF (ci > 1) ci = 1
253 cilatenight = ci
254
255 !Estimation of radD and radI if not measured after Reindl et al. (1990)
256 IF (onlyglobal == 1) THEN
257 CALL diffusefraction(radg, altitude, kt, ta, rh, radi, radd)
258 END IF
259
260 !Shadow images
261 IF (usevegdem == 1) THEN ! with vegetation
262 ! CALL shadowingfunction_20(azimuth, altitude, scale, amaxvalue)
263 ! shadow = (sh - (1 - vegsh)*(1 - psi))
264 ELSE ! without vegetation
265 CALL shadowingfunction_urban(azimuth, altitude, scale, shadow)
266 vegsh = 1.0d0
267 ! shadow = sh
268 END IF
269
270 !Ground View Factors based on shadow patterns and sunlit walls
271 gvf = 0.0d0
272 ! CALL wallinsun_veg(azimuth,sunwall)
273 DO j = 1, SIZE(azimutha)
274 CALL sunonsurface_veg(azimutha(j), scale, buildings, first, second, psi, sos)
275 gvf = gvf + sos
276 END DO
277 gvf = gvf/SIZE(azimutha) + (buildings*(-1) + 1)
278
279 !Building height angle from svf
280 CALL cylindric_wedge(zen, svfalfa, f_sh) !Fraction shadow on building walls based on sun altitude and svf
281 !F_sh(isnan(F_sh))=0.5 FIXTHIS
282
283!!! Calculation of shortwave daytime radiative fluxes !!!
284 kdown2d = radi*shadow*sin(altitude*deg2rad) &
285 + radd*svf_blgd_veg &
286 + radg*alb_bldg*(1 - svf_blgd_veg)*(1 - f_sh) !*sin(altitude(i)*(pi/180));
287
288 kup2d = alb_ground*( &
289 radi*gvf*sin(altitude*deg2rad) &
290 + radd*svf_blgd_veg &
291 + radg*alb_bldg*(1 - svf_blgd_veg)*(1 - f_sh))
292
293 CALL kside_veg_v24( &
294 shadow, f_sh, &
295 radi, radg, radd, azimuth, altitude, psi, t, alb_ground, & ! input
296 keast, knorth, ksouth, kwest) ! output
297
298!!!! Surface temperature parameterisation during daytime !!!!
299 ! dfm = ABS(172 - DOY) !Day from midsommer
300 ! Tgamp = 0.000006*dfm**3 - 0.0017*dfm**2 + 0.0127*dfm + 17.084 + Tstart !sinus function for daily surface temperature wave
301 ! !Tg=Tgamp*sin(((hour-rise)/(15-rise))*pi/2)-Tstart ! check if this should be 15 (3 pm)
302 ! Tg = Tgamp*SIN(((dectime - DOY - SNUP/24)/(15./24 - SNUP/24))*pi/2) - Tstart !new sunrise time 2014a
303 ! IF (Tg < 0) Tg = 0! temporary for removing low Tg during morning 20140513, from SOLWEIG1D
304 ! s = (dectime - DOY - SNUP/24)
305 ! !New estimation of Tg reduction for non-clear situation based on Reindl et al. 1990
306 ! Ktc = 1.0
307 ! CALL diffusefraction(I0, altitude, Ktc, Ta, RH, radI0, s)
308 ! corr = 0.1473*LOG(90.-(zen*RAD2DEG)) + 0.3454 ! 20070329 temporary correction of latitude from Lindberg et al. 2008
309 ! CI_Tg = (radI/radI0) + (1.-corr)
310
311 ! IF (CI_Tg > 1) CI_Tg = 1 !!FIX THIS?? .and. CI_Tg<inf then CI_Tg=1
312
313 ! Tg = Tg*CI_Tg !new estimation
314 tw = tg
315
316!!!! Lup, daytime !!!!
317 !Surface temperature wave delay - new as from 2014a
318 tgmap0 = gvf*tg + ta ! current timestep
319 tgmap1 = tgmap0 !"first in morning"
320 ! IF (firstdaytime == 1) THEN !"first in morning"
321 ! Tgmap1 = Tgmap0
322 ! END IF
323 ! IF (timeadd >= (59/1440.)) THEN !more or equal to 59 min
324 ! ! weight1 = EXP(-33.27*timeadd) !surface temperature delay function - 1 step
325 ! ! Tgmap1 = Tgmap0*(1 - weight1) + Tgmap1*weight1
326 ! ! Lup2d = SBC*emis_ground*((Tgmap1 + 273.15)**4)
327 ! timestepdec=tstep/(60*60*24)
328 ! IF (timestepdec > (59/1440.)) THEN
329 ! timeadd = timestepdec
330 ! ELSE
331 ! timeadd = 0
332 ! END IF
333 ! ELSE
334 ! timeadd = timeadd + timestepdec
335 ! END IF
336 timeadd = (dectime + 1) - doy - (it/24.)
337 weight1 = exp(-33.27*timeadd) !surface temperature delay function - 1 step
338 tgmap1 = tgmap0*(1 - weight1) + tgmap1*weight1
339 lup2d = sbc*emis_ground*((tgmap1 + 273.15)**4)
340 ! firstdaytime = 0
341
342 ELSE !!!!!!! NIGHTTIME !!!!!!!!
343
344 !Nocturnal cloud fraction from Offerle et al. 2003
345 IF (dectime < (doy + 0.5) .AND. dectime > doy .AND. altitude < 1.0) THEN !! THIS NEED SOME THOUGHT 20150211
346 !j=0
347 !do while (dectime<(DOY+SNUP/24))
348 !! call ConvertMetData(ith+j) ! read data at sunrise ??
349 ! j=j+1
350 !end do
351 !call NARP_cal_SunPosition(year,dectime,timezone,lat,lng,alt,azimuth,zenith_deg)!this is not good
352 !zen=zenith_deg*DEG2RAD
353 !call clearnessindex_2013b(zen,DOY,Temp_C,RH/100,avkdn,lat,Press_hPa,I0,CI,Kt,I0et,CIuncorr)
354 !!call ConvertMetData(ith) ! read data at current timestep again ??
355 !call NARP_cal_SunPosition(year,dectime,timezone,lat,lng,alt,azimuth,zenith_deg)!this is not good
356
357 ci = 1.0
358 ELSE
359 ! IF (SolweigCount == 1) THEN
360 ! CI = 1.0
361 ! ELSE
362 ! CI = CIlatenight
363 ! END IF
364 ci = cilatenight
365 END IF
366
367 tw = 0.0
368 ! Tg = 0.0
369 radg = 0
370 radi = 0
371 radd = 0
372
373 !Nocturnal Kfluxes set to 0
374 kdown2d = 0.0
375 kwest = 0.0
376 kup2d = 0.0
377 keast = 0.0
378 ksouth = 0.0
379 knorth = 0.0
380 shadow = 0.0
381 gvf = 0.0
382
383!!! Lup !!!
384 lup2d = sbc*emis_ground*((knight + ta + tg + 273.15)**4)
385 firstdaytime = 1
386 END IF
387
388!!! Ldown !!!
389 IF (solweig_ldown == 1) THEN ! fixed for non-clear situations 20140701
390 ldown2d = (svf + svfveg - 1)*emis_sky*sbc*((ta + 273.15)**4) &
391 + (2 - svfveg - svfaveg)*emis_wall*sbc*((ta + 273.15)**4) &
392 + (svfaveg - svf)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
393 + (2 - svf - svfveg)*(1 - emis_wall)*emis_sky*sbc*((ta + 273.15)**4) !Jonsson et al. (2006)
394 !Ldown2d=Ldown2d-25 ! Shown by Jonsson et al. (2006) and Duarte et al. (2006) ! Removed as from 20140701
395 IF (ci > 1) ci = 1 !!FIX THIS?? .and. CI<inf) CI=1
396 !if (CI < 0.95) then !non-clear conditions
397 c = 1 - ci
398 !Ldown2d=Ldown2d*(1-c)+c*SBC*((Ta+273.15)**4)
399 ldown2d = ldown2d*(1 - c) + c*( &
400 (svf + svfveg - 1)*sbc*((ta + 273.15)**4) &
401 + (2 - svfveg - svfaveg)*emis_wall*sbc*((ta + 273.15)**4) &
402 + (svfaveg - svf)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
403 + (2 - svf - svfveg)*(1 - emis_wall)*sbc*((ta + 273.15)**4))
404 !end if
405 ELSE
406 ldown2d = (svf + svfveg - 1)*ldown &
407 + (2 - svfveg - svfaveg)*emis_wall*sbc*((ta + 273.15)**4) &
408 + (svfaveg - svf)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
409 + (2 - svf - svfveg)*(1 - emis_wall)*ldown !Jonsson et al. (2006)
410 END IF
411
412!!! Lside !!!
413 CALL lside_veg_v2(ldown2d, lup2d, &
414 altitude, ta, tw, sbc, emis_wall, emis_sky, t, ci, azimuth, zen, ldown, svfalfa, &
415 least, lnorth, lsouth, lwest)
416
417!!! Calculation of radiant flux density and Tmrt !!!
418 sstr = absk*(kdown2d*fup + kup2d*fup + knorth*fside + keast*fside + ksouth*fside + kwest*fside) &
419 + absl*(ldown2d*fup + lup2d*fup + lnorth*fside + least*fside + lsouth*fside + lwest*fside)
420 tmrt = sqrt(sqrt((sstr/(absl*sbc)))) - 273.15
421
422 ! IF (SOLWEIGpoi_out == 1) THEN
423 ! dataOutSOLWEIG(SolweigCount, 1:4, iMBi) = [iy, id, it, imin]
424 ! dataOutSOLWEIG(SolweigCount, 5, iMBi) = dectime
425 ! dataOutSOLWEIG(SolweigCount, 6:ncolumnsdataOutSOL, iMBi) = &
426 ! [azimuth, altitude, radG, radD, radI, &
427 ! Kdown2d(row, col), Kup2d(row, col), Ksouth(row, col), Kwest(row, col), Knorth(row, col), Keast(row, col), &
428 ! Ldown2d(row, col), Lup2d(row, col), Lsouth(row, col), Lwest(row, col), Lnorth(row, col), Least(row, col), &
429 ! Tmrt(row, col), I0, CI, gvf(row, col), shadow(row, col), svf(row, col), svf_blgd_veg(row, col), Ta, Tg]
430 ! END IF
431 dataoutlinesolweig = [azimuth, altitude, radg, radi, radd, &
432 kdown2d(row, col), kup2d(row, col), &
433 ksouth(row, col), kwest(row, col), knorth(row, col), keast(row, col), &
434 ldown2d(row, col), lup2d(row, col), &
435 lsouth(row, col), lwest(row, col), lnorth(row, col), least(row, col), &
436 tmrt(row, col), i0, ci, gvf(row, col), shadow(row, col), svf(row, col), &
437 svf_blgd_veg(row, col), ta, tg]
438
439 ! CALL SaveGrids
440
441 DEALLOCATE (tmp)
442 DEALLOCATE (knight)
443 DEALLOCATE (tgmap0)
444 DEALLOCATE (svf_blgd_veg)
445
446 ! external grids
447 ! DEALLOCATE (Kdown2d)
448 ! DEALLOCATE (Kup2d)
449 ! DEALLOCATE (Knorth)
450 ! DEALLOCATE (Kwest)
451 ! DEALLOCATE (Ksouth)
452 ! DEALLOCATE (Keast)
453 ! DEALLOCATE (Ldown2d)
454 ! DEALLOCATE (Lup2d)
455 ! DEALLOCATE (Lnorth)
456 ! DEALLOCATE (Lwest)
457 ! DEALLOCATE (Lsouth)
458 ! DEALLOCATE (Least)
459
460 END SUBROUTINE solweig_cal_main
461
462 FUNCTION cal_ratio_height2width(lamdaP, lamdaF) RESULT(HW)
463 IMPLICIT NONE
464 REAL(kind(1d0)), PARAMETER :: a = 0.5598
465 REAL(kind(1d0)), PARAMETER :: b = -0.2485
466 REAL(kind(1d0)), PARAMETER :: c = 0.4112
467 REAL(kind(1d0)), PARAMETER :: d = -0.02528
468
469 REAL(kind(1d0)), INTENT(in) :: lamdap ! plan area fraction
470 REAL(kind(1d0)), INTENT(in) :: lamdaf ! frontal area fraction
471 REAL(kind(1d0)) :: hw
472 REAL(kind(1d0)) :: lamdaw !wall area fraction (wallarea / total area)
473
474 lamdaw = 4*lamdaf ! assuming square shaped buildings
475
476 hw = (lamdaw*lamdap)/(2*lamdap*(1 - lamdap))
477
478 END FUNCTION cal_ratio_height2width
479
480 FUNCTION hwtosvf_ground(hw) RESULT(svfGround)
481 IMPLICIT NONE
482 REAL(kind(1d0)), PARAMETER :: a = 0.5598
483 REAL(kind(1d0)), PARAMETER :: b = -0.2485
484 REAL(kind(1d0)), PARAMETER :: c = 0.4112
485 REAL(kind(1d0)), PARAMETER :: d = -0.02528
486
487 REAL(kind(1d0)), INTENT(in) :: hw
488 REAL(kind(1d0)) :: svfground
489
490 ! SvfGround: Parameterisation based on NYC data (500x500 meter grid)
491 svfground = a*exp(b*hw) + c*exp(d*hw)
492
493 END FUNCTION hwtosvf_ground
494
495 FUNCTION hwtosvf_roof(hw) RESULT(svfRoof)
496 IMPLICIT NONE
497 REAL(kind(1d0)), PARAMETER :: e = 0.5572
498 REAL(kind(1d0)), PARAMETER :: f = 0.0589
499 REAL(kind(1d0)), PARAMETER :: g = 0.4143
500
501 REAL(kind(1d0)), INTENT(in) :: hw
502 REAL(kind(1d0)) :: svfroof
503
504 ! SvfGround: Parameterisation based on NYC data (500x500 meter grid)
505 svfroof = e*exp(-f*hw) + g
506
507 END FUNCTION hwtosvf_roof
508
510 SUBROUTINE clearnessindex_2013b(zen, DOY, Ta, RH, radG, lat, P_kPa, I0, CI, Kt, I0et, CIuncorr)
511 !Last modified:
512 !LJ 27 Jan 2016 - Removal of tabs
513
514 IMPLICIT NONE
515 ! Use somemodule
516 INTEGER, INTENT(in) :: DOY
517
518 REAL(KIND(1D0)), INTENT(in) :: zen
519 REAL(KIND(1D0)), INTENT(in) :: Ta
520 REAL(KIND(1D0)), INTENT(in) :: RH
521 REAL(KIND(1D0)), INTENT(in) :: P_kPa
522 REAL(KIND(1D0)), INTENT(in) :: radG
523 REAL(KIND(1D0)), INTENT(in) :: lat
524 REAL(KIND(1D0)), INTENT(out) :: I0et
525 REAL(KIND(1D0)), INTENT(out) :: CIuncorr
526 REAL(KIND(1D0)), INTENT(out) :: CI
527 REAL(KIND(1D0)), INTENT(out) :: I0
528 REAL(KIND(1D0)), INTENT(out) :: Kt
529 REAL(KIND(1D0)) :: iG, Itoa, p
530 REAL(KIND(1D0)), DIMENSION(4) :: G
531 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
532
533 ! Variable declarations
534 real*8 :: a2
535 !logical :: b !>
536 real*8 :: b2
537 real*8 :: corr
538 real*8 :: d
539 real*8 :: m
540 real*8 :: tar
541 real*8 :: td
542 real*8 :: trpg
543 real*8 :: tw
544 real*8 :: u
545 !
546
547 ! Clearness Index at the Earth's surface calculated from Crawford and Duchon 1999
548
549 IF (p_kpa == -999) THEN
550 p = 1013 !Pressure in millibars
551 ELSE
552 p = p_kpa*10 !Convert from hPa to millibars
553 END IF
554 !Effective solar constant
555 itoa = 1370
556 !call solar_ESdist(jday,D)
557 CALL sun_distance(doy, d)
558 !D=sun_distance(jday) !irradiance differences due to Sun-Earth distances
559 m = 35.*cos(zen)*((1224.*(cos(zen)**2) + 1.)**(-1./2.)) !optical air mass at p=1013
560 trpg = 1.021 - 0.084*(m*(0.000949*p + 0.051))**0.5 !Transmission coefficient for Rayliegh scattering and permanent gases
561
562 ! empirical constant depending on latitude
563 g = 0
564 IF (lat < 10) THEN
565 g = [3.37, 2.85, 2.80, 2.64]
566 ELSE IF (lat >= 10 .AND. lat < 20) THEN
567 g = [2.99, 3.02, 2.70, 2.93]
568 ELSE IF (lat >= 20 .AND. lat < 30) THEN
569 g = [3.60, 3.00, 2.98, 2.93]
570 ELSE IF (lat >= 30 .AND. lat < 40) THEN
571 g = [3.04, 3.11, 2.92, 2.94]
572 ELSE IF (lat >= 40 .AND. lat < 50) THEN
573 g = [2.70, 2.95, 2.77, 2.71]
574 ELSE IF (lat >= 50 .AND. lat < 60) THEN
575 g = [2.52, 3.07, 2.67, 2.93]
576 ELSE IF (lat >= 60 .AND. lat < 70) THEN
577 g = [1.76, 2.69, 2.61, 2.61]
578 ELSE IF (lat >= 70 .AND. lat < 80) THEN
579 g = [1.60, 1.67, 2.24, 2.63]
580 ELSE IF (lat >= 80 .AND. lat < 90) THEN
581 g = [1.11, 1.44, 1.94, 2.02]
582 END IF
583 IF (doy > 335 .OR. doy <= 60) THEN
584 ig = g(1)
585 ELSE IF (doy > 60 .AND. doy <= 152) THEN
586 ig = g(2)
587 ELSE IF (doy > 152 .AND. doy <= 244) THEN
588 ig = g(3)
589 ELSE IF (doy > 244 .AND. doy <= 335) THEN
590 ig = g(4)
591 END IF
592 !dewpoint calculation
593 a2 = 17.27
594 b2 = 237.7
595 td = (b2*(((a2*ta)/(b2 + ta)) + log(rh)))/(a2 - (((a2*ta)/(b2 + ta)) + log(rh)))
596 td = (td*1.8) + 32 !Dewpoint (degF)
597 u = exp(0.1133 - log(ig + 1) + 0.0393*td) !Precipitable water
598 tw = 1 - 0.077*((u*m)**0.3) !Transmission coefficient for water vapor
599 tar = 0.935**m !Transmission coefficient for aerosols
600
601 i0 = itoa*cos(zen)*trpg*tw*d*tar
602
603!!! This needs to be checked !!!
604 !b=I0==abs(zen)>pi/2
605 !I0(b==1)=0
606 !clear b
607 !if (not(isreal(I0))) then
608 ! I0=0
609 !end if
610
611 corr = 0.1473*log(90 - (zen/pi*180)) + 0.3454 ! 20070329
612
613 ciuncorr = radg/i0
614 ci = ciuncorr + (1 - corr)
615 i0et = itoa*cos(zen)*d !extra terrestial solar radiation
616 kt = radg/i0et
617
618!!!! This needs to be checked !!!
619 !if (isnan(CI)) then
620 ! CI=NaN
621 !end if
622 END SUBROUTINE clearnessindex_2013b
623
624 !===============================================================================
625
626 SUBROUTINE sun_distance(jday, D)
627
628 ! Calculates solar irradiance variation based on mean earth sun distance
629 ! with day of year as input.
630 ! Partridge and Platt, 1975
631
632 INTEGER :: jday
633 REAL(KIND(1D0)) :: b, D
634
635 b = 2*3.141592654*jday/365
636 d = sqrt(1.00011 + 0.034221*cos(b) + 0.001280*sin(b) + 0.000719*cos(2*b) + 0.000077*sin(2*b))
637
638 END SUBROUTINE sun_distance
639
640 SUBROUTINE cylindric_wedge(zen, svfalfa, F_sh)
641 ! Fraction of sunlit walls based on sun altitude and svf wieghted building angles
642
643 IMPLICIT NONE
644
645 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
646 REAL(KIND(1D0)), INTENT(in) :: zen
647 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: svfalfa
648 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: F_sh
649
650 REAL(KIND(1D0)) :: beta
651 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: alfa, xa, ha, hkil, ba
652 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: Ai, phi, qa, Za
653 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: ukil, Ssurf
654 !real(kind(1d0)), dimension(1,1) ::
655 !real(kind(1d0)), dimension(1,1) ::
656 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: sos, Tgmap1
657 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: gvf, Tmrt, shadow, Sstr, sunwall
658
659 ALLOCATE (alfa(1, 1))
660 ALLOCATE (ba(1, 1))
661 ALLOCATE (ha(1, 1))
662 ALLOCATE (xa(1, 1))
663 ALLOCATE (qa(1, 1))
664 ALLOCATE (za(1, 1))
665 ALLOCATE (phi(1, 1))
666 ALLOCATE (ukil(1, 1))
667 ALLOCATE (ai(1, 1))
668 ALLOCATE (ssurf(1, 1))
669 ALLOCATE (hkil(1, 1))
670
671 beta = zen
672 alfa = svfalfa
673
674 xa = 1.-2./(tan(alfa)*tan(beta))
675 ha = 2./(tan(alfa)*tan(beta))
676 ba = (1./tan(alfa))
677 hkil = 2.*ba*ha
678
679 qa = 0.0d0
680
681 WHERE (xa < 0) !qa(xa<0)=tan(beta)/2
682 qa = tan(beta)/2
683 END WHERE
684
685 za = 0.0d0
686
687 phi = 0.0d0
688
689 ai = 0.0d0
690
691 ukil = 0.0d0
692 WHERE (xa < 0)
693 !Za(xa<0)=((ba(xa<0).**2)-((qa(xa<0).**2)./4)).**0.5
694 za = (ba**2 - qa**2/4.)**0.5
695 !phi(xa<0)=atan(Za(xa<0)./qa(xa<0))
696 phi = atan(za/qa)
697 !A1(xa<0)=(sin(phi(xa<0))-phi(xa<0).*cos(phi(xa<0)))./(1-cos(phi(xa<0)))
698 ai = (sin(phi) - phi*cos(phi))/(1 - cos(phi))
699 !ukil(xa<0)=2*ba(xa<0).*xa(xa<0).*A1(xa<0)
700 ukil = 2*ba*xa*ai
701 END WHERE
702
703 ssurf = hkil + ukil
704
705 f_sh = (2*pi*ba - ssurf)/(2*pi*ba) !Xa
706
707 DEALLOCATE (alfa)
708 DEALLOCATE (ba)
709 DEALLOCATE (ha)
710 DEALLOCATE (xa)
711 DEALLOCATE (qa)
712 DEALLOCATE (za)
713 DEALLOCATE (phi)
714 DEALLOCATE (ukil)
715 DEALLOCATE (ai)
716 DEALLOCATE (ssurf)
717 DEALLOCATE (hkil)
718
719 END SUBROUTINE cylindric_wedge
720
721 ! This subroutine estimates diffuse and directbeam radiation according to
722 ! Reindl et al (1990), Solar Energy 45:1
723 SUBROUTINE diffusefraction(radG, altitude, Kt, Ta, RH, radI, radD)
724 IMPLICIT NONE
725
726 REAL(KIND(1D0)), INTENT(in) :: radG
727 REAL(KIND(1D0)), INTENT(in) :: altitude
728 REAL(KIND(1D0)), INTENT(in) :: Kt
729 REAL(KIND(1D0)), INTENT(in) :: Ta
730 REAL(KIND(1D0)), INTENT(in) :: RH
731 REAL(KIND(1D0)), INTENT(out) :: radD ! direct radiation
732 REAL(KIND(1D0)), INTENT(out) :: radI ! diffusive radiation
733 REAL(KIND(1D0)) :: alfa
734 ! REAL(KIND(1D0)), PARAMETER :: DEG2RAD = 0.017453292, RAD2DEG = 57.29577951 !!Already defined in AllocateArray module. Delete??
735
736 alfa = altitude*deg2rad
737
738 IF (ta <= -99 .OR. rh <= -99) THEN !.or. isnan(Ta) .or. isnan(RH)) then
739 IF (kt <= 0.3) THEN
740 radd = radg*(1.020 - 0.248*kt)
741 ELSE IF (kt > 0.3 .AND. kt < 0.78) THEN
742 radd = radg*(1.45 - 1.67*kt)
743 ELSE IF (kt >= 0.78) THEN
744 radd = radg*0.147
745 END IF
746 ELSE
747 !RH=RH/100
748 IF (kt <= 0.3) THEN
749 radd = radg*(1 - 0.232*kt + 0.0239*sin(alfa) - 0.000682*ta + 0.0195*(rh/100))
750 ELSE IF (kt > 0.3 .AND. kt < 0.78) THEN
751 radd = radg*(1.329 - 1.716*kt + 0.267*sin(alfa) - 0.00357*ta + 0.106*(rh/100))
752 ELSE IF (kt >= 0.78) THEN
753 radd = radg*(0.426*kt - 0.256*sin(alfa) + 0.00349*ta + 0.0734*(rh/100))
754 END IF
755 END IF
756 ! correction of radD
757 radd = max(0.d0, radd)
758 radd = min(radg, radd)
759
760 ! calculation of diffuse radiation
761 radi = (radg - radd)/(sin(alfa))
762
763 !! Corrections for low sun altitudes (20130307)
764 ! IF (radI < 0) THEN
765 ! radI = 0
766 ! END IF
767 ! IF (radD < 0) THEN
768 ! radD = 0
769 ! END IF
770
771 IF (altitude < 1 .AND. radi > radg) THEN
772 radi = radg
773 END IF
774
775 ! IF (radD > radG) THEN
776 ! radD = radG
777 ! END IF
778 END SUBROUTINE diffusefraction
779 ! This subroutine loads a ESRIACSII grid as a 2D array
780 ! Last modified:
781 ! LJ 27 Jan 2016-removal of tabs
782 !-----------------------------------------------------
783 ! SUBROUTINE LoadEsriAsciiGrid(GridPath, GridName, xllcornerlo, yllcornerlo, cellsizelo, NoDatalo)
784
785 ! IMPLICIT NONE
786 ! REAL(KIND(1d0)) :: xllcornerlo, yllcornerlo, cellsizelo, NoDatalo
787 ! INTEGER :: col, row
788 ! CHARACTER(len=100) :: GridPath, GridName, GridFile, n
789 ! !real(kind(1d0)),intent(out) :: tempgrid
790 ! !real(kind(1d0)),dimension(1,1),intent(out) :: tempgrid!,allocatable
791
792 ! ! Loading DSM
793 ! !GridPath='D:\SOLWEIG2013b_Fortran\Inputdata\'
794 ! !GridName='kr_dem.asc'
795 ! GridFile = TRIM(GridPath)//TRIM(GridName)
796 ! OPEN (99, File=GridFile, status='old')
797
798 ! ! Read Header
799 ! READ (99, *) n, 1
800 ! READ (99, *) n, 1
801 ! READ (99, *) n, xllcornerlo
802 ! READ (99, *) n, yllcornerlo
803 ! READ (99, *) n, cellsizelo
804 ! READ (99, *) n, NoDatalo
805
806 ! ALLOCATE (tempgrid(1, 1))
807
808 ! ! Read Matrix
809 ! DO row = 1, 1
810 ! READ (99, *) (tempgrid(row, col), col=1, 1)
811 ! END DO
812 ! CLOSE (99)
813
814 ! RETURN
815 ! END SUBROUTINE LoadEsriAsciiGrid
816
817! ! This subroutine saves a 2D array as an ESRIACSII grid
818! SUBROUTINE SaveEsriAsciiGrid(GridPath, GridName, xllcornerlo, yllcornerlo, cellsizelo, NoDatalo)
819
820! IMPLICIT NONE
821! REAL(KIND(1d0)) :: xllcornerlo, yllcornerlo, cellsizelo, NoDatalo
822! INTEGER :: col, row
823! CHARACTER(len=100) :: GridPath, GridName, GridFile
824! !integer :: 1,1!,intent(in)
825! !real(kind(1d0)), allocatable, dimension(:,:):: grid
826! ! real(kind(1d0)),dimension(1,1) :: grid!,allocatable
827
828! ! Loading DSM
829! !GridPath='D:\SOLWEIG2013b_Fortran\Inputdata\'
830! !GridName='kr_dem.asc'
831! GridFile = TRIM(GridPath)//TRIM(GridName)
832! OPEN (94, File=GridFile, status='unknown')
833
834! ! Read Header
835! WRITE (94, "(A5,1x,I0)") 'ncols', 1
836! WRITE (94, "(A5,1x,I0)") 'nrows', 1
837! WRITE (94, "(A9,1x,F0.2)") 'xllcorner', xllcornerlo
838! WRITE (94, "(A9,1x,F0.2)") 'yllcorner', yllcornerlo
839! WRITE (94, "(A8,1x,F0.2)") 'cellsize', cellsizelo
840! WRITE (94, "(A12,1x,F0.2)") 'NODATA_value', NoDatalo
841
842! ! write Matrix
843! DO row = 1, 1
844! WRITE (94, 100) (savegrid(row, col), col=1, 1)
845! END DO
846! CLOSE (94)
847! 100 FORMAT(200(f6.2, 1x))
848
849! RETURN
850! END SUBROUTINE SaveEsriAsciiGrid
851
852 ! setup for SOLWEIG
853 ! FL may 2014
854! SUBROUTINE SOLWEIG_Initial
855! !Last modified LJ 27 Jan 2016 - Removal of Tabs and real-int fixes
856
857! ! USE matsize ! All allocatable grids and related variables used in SOLWEIG
858
859! IMPLICIT NONE
860
861! CHARACTER(len=100) :: Path, GridFile, GridFolder
862! REAL(KIND(1d0)) :: vegmax
863! CHARACTER(len=100), DIMENSION(5) :: svfname
864! CHARACTER(len=100), DIMENSION(10):: svfvegname
865! LOGICAL :: exist
866! INTEGER :: firstday
867
868! NAMELIST /SOLWEIGinput/ Posture, & ! 1.Standing, 2.Sitting
869! absL, & ! Absorption coefficient of longwave radiation of a person
870! absK, & ! Absorption coefficient of shortwave radiation of a person
871! heightgravity, & ! Center of gravity for a standing person
872! usevegdem, & ! With vegetation (1)
873! DSMPath, & ! Path to DSMs
874! DSMname, & ! Ground and building DSM
875! CDSMname, & ! Canopy DSM
876! TDSMname, & ! Trunk zone DSM
877! TransMin, & ! Transmissivity of K through decidious vegetation (leaf on)
878! TransMax, & ! Transmissivity of K through decidious vegetation (leaf off)
879! SVFPath, & ! Path to SVFs
880! SVFsuffix, & !
881! buildingsname, & ! Boolean matrix for locations of building pixels
882! row, & ! X coordinate for point of interest
883! col, & ! Y coordinate for point of interest
884! onlyglobal, & ! if no diffuse and direct, then =1
885! SOLWEIGpoi_out, & ! write output variables at point of interest
886! Tmrt_out, & ! write Tmrt grid to file
887! Lup2d_out, & ! write Lup grid to file
888! Ldown2d_out, & ! write Ldown grid to file
889! Kup2d_out, & ! write Kup grid to file
890! Kdown2d_out, & ! write Kdown grid to file
891! GVF_out, & ! write GroundViewFactor grid to file
892! SOLWEIG_ldown, & ! 1= use SOLWEIG code to estimate Ldown, 0=use SEUWS
893! OutInterval, & ! Output interval in minutes
894! RunForGrid ! If only one grid should be run. All grids -999
895
896! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
897! !Read in the SOLWEIGinput.nml file
898! OPEN (52, file=TRIM(FileInputPath)//'SOLWEIGinput.nml', err=274, status='old')
899! READ (52, nml=SOLWEIGinput)
900! CLOSE (52)
901
902! SolweigCount = 1
903
904! IF (OutInterval == 60) THEN
905! OutInterval = 0
906! ENDIF
907
908! IF (Posture == 1) THEN
909! Fside = 0.22
910! Fup = 0.06
911! ELSE
912! Fside = 0.1666666667
913! Fup = 0.166666667
914! ENDIF
915
916! !timestepdec=real(t_interval)/(3600*24)
917! timestepdec = REAL(OutInterval)/(REAL(t_interval)*24.)
918
919! !!! Loading DSM !!!
920! Path = TRIM(FileInputPath)//TRIM(DSMPath)
921! CALL LoadEsriAsciiGrid(Path, DSMName, xllcorner, yllcorner, cellsize, NoData)
922! ALLOCATE (a(1, 1))
923! a = tempgrid
924! DEALLOCATE (tempgrid)
925
926! scale = 1/cellsize
927
928! GridFolder = TRIM(FileOutputPath)//'Grids'
929
930! ! Create grid folder !! This does not work in windows. needs to be done in python
931! !inquire(file=GridFolder, exist=exist)
932! !if (exist) then
933! !else
934! ! makedirectory = 'mkdir ' // trim(GridFolder)
935! ! call system(makedirectory)
936! !end if
937
938! !!! Set up for vegetation scheme, or not !!!
939! IF (usevegdem == 1) THEN
940! ! Calculating transmissivity of short wave radiation through vegetation based on decid LAI
941! transperLAI = (TransMax - TransMin)/(LAImax(2) - LAImin(2))
942! firstday = INT(MetForcingData(1, 2, 1))
943! trans = TransMin + (LAImax(2) - LAI(firstday - 1, 2))*transperLAI
944
945! ! Loading vegDSM (SDSM)
946! Path = TRIM(FileInputPath)//TRIM(DSMPath)
947! CALL LoadEsriAsciiGrid(Path, CDSMname, xllcorner, yllcorner, cellsize, NoData)
948! ALLOCATE (vegdem(1, 1))
949! vegdem = tempgrid
950! DEALLOCATE (tempgrid)
951
952! ! Loading trunkDSM (TDSM)
953! Path = TRIM(FileInputPath)//TRIM(DSMPath)
954! CALL LoadEsriAsciiGrid(Path, TDSMname, xllcorner, yllcorner, cellsize, NoData)
955! ALLOCATE (vegdem2(1, 1))
956! vegdem2 = tempgrid
957! DEALLOCATE (tempgrid)
958
959! ! amaxvalue (used in calculation of vegetation shadows)
960! vegmax = MAXVAL(vegdem)
961! amaxvalue = MAXVAL(a) - MINVAL(a)
962! amaxvalue = MAX(amaxvalue, vegmax)
963
964! ! Elevation vegdems if buildingDSM includes ground heights
965! vegdem = vegdem + a
966! WHERE (vegdem == a)
967! vegdem = 0.0
968! END WHERE
969! vegdem2 = vegdem2 + a;
970! WHERE (vegdem2 == a)
971! vegdem2 = 0.0
972! END WHERE
973! ! Bush separation
974! ALLOCATE (bush(1, 1))
975! WHERE ((vegdem > 0) .AND. (vegdem2 == 0))
976! bush = vegdem
977! ELSEWHERE
978! bush = 0.0
979! END WHERE
980! ELSE
981! trans = 1.00;
982! ENDIF
983
984! !!! Loading/creating SVFs !!!
985! Path = TRIM(FileInputPath)//TRIM(SVFPath)//TRIM(SVFsuffix)
986! svfname = ['svf.asc ', 'svfE.asc', 'svfN.asc', 'svfW.asc', 'svfS.asc']
987! svfvegname = ['svfveg.asc ', 'svfEveg.asc ', 'svfNveg.asc ', 'svfWveg.asc ', 'svfSveg.asc ', &
988! 'svfaveg.asc ', 'svfEaveg.asc', 'svfNaveg.asc', 'svfWaveg.asc', 'svfSaveg.asc']
989! ! SVFs, Should be done as a loop... ! How to change variable in a loop???
990! CALL LoadEsriAsciiGrid(Path, svfname(1), xllcorner, yllcorner, cellsize, NoData)
991! ALLOCATE (svf(1, 1)); svf = tempgrid; DEALLOCATE (tempgrid)
992! CALL LoadEsriAsciiGrid(Path, svfname(2), xllcorner, yllcorner, cellsize, NoData)
993! ALLOCATE (svfE(1, 1)); svfE = tempgrid; DEALLOCATE (tempgrid)
994! CALL LoadEsriAsciiGrid(Path, svfname(3), xllcorner, yllcorner, cellsize, NoData)
995! ALLOCATE (svfN(1, 1)); svfN = tempgrid; DEALLOCATE (tempgrid)
996! CALL LoadEsriAsciiGrid(Path, svfname(4), xllcorner, yllcorner, cellsize, NoData)
997! ALLOCATE (svfW(1, 1)); svfW = tempgrid; DEALLOCATE (tempgrid)
998! CALL LoadEsriAsciiGrid(Path, svfname(5), xllcorner, yllcorner, cellsize, NoData)
999! ALLOCATE (svfS(1, 1)); svfS = tempgrid; DEALLOCATE (tempgrid)
1000! CALL LoadEsriAsciiGrid(Path, svfvegname(1), xllcorner, yllcorner, cellsize, NoData)
1001! ALLOCATE (svfveg(1, 1)); svfveg = tempgrid; DEALLOCATE (tempgrid)
1002! CALL LoadEsriAsciiGrid(Path, svfvegname(2), xllcorner, yllcorner, cellsize, NoData)
1003! ALLOCATE (svfEveg(1, 1)); svfEveg = tempgrid; DEALLOCATE (tempgrid)
1004! CALL LoadEsriAsciiGrid(Path, svfvegname(3), xllcorner, yllcorner, cellsize, NoData)
1005! ALLOCATE (svfNveg(1, 1)); svfNveg = tempgrid; DEALLOCATE (tempgrid)
1006! CALL LoadEsriAsciiGrid(Path, svfvegname(4), xllcorner, yllcorner, cellsize, NoData)
1007! ALLOCATE (svfWveg(1, 1)); svfWveg = tempgrid; DEALLOCATE (tempgrid)
1008! CALL LoadEsriAsciiGrid(Path, svfvegname(5), xllcorner, yllcorner, cellsize, NoData)
1009! ALLOCATE (svfSveg(1, 1)); svfSveg = tempgrid; DEALLOCATE (tempgrid)
1010! CALL LoadEsriAsciiGrid(Path, svfvegname(6), xllcorner, yllcorner, cellsize, NoData)
1011! ALLOCATE (svfaveg(1, 1)); svfaveg = tempgrid; DEALLOCATE (tempgrid)
1012! CALL LoadEsriAsciiGrid(Path, svfvegname(7), xllcorner, yllcorner, cellsize, NoData)
1013! ALLOCATE (svfEaveg(1, 1)); svfEaveg = tempgrid; DEALLOCATE (tempgrid)
1014! CALL LoadEsriAsciiGrid(Path, svfvegname(8), xllcorner, yllcorner, cellsize, NoData)
1015! ALLOCATE (svfNaveg(1, 1)); svfNaveg = tempgrid; DEALLOCATE (tempgrid)
1016! CALL LoadEsriAsciiGrid(Path, svfvegname(9), xllcorner, yllcorner, cellsize, NoData)
1017! ALLOCATE (svfWaveg(1, 1)); svfWaveg = tempgrid; DEALLOCATE (tempgrid)
1018! CALL LoadEsriAsciiGrid(Path, svfvegname(10), xllcorner, yllcorner, cellsize, NoData)
1019! ALLOCATE (svfSaveg(1, 1)); svfSaveg = tempgrid; DEALLOCATE (tempgrid)
1020
1021! !!! Loading buildings grid !!!
1022! Path = TRIM(FileInputPath)//TRIM(DSMPath)
1023! GridFile = TRIM(Path)//TRIM(buildingsname)
1024! INQUIRE (file=GridFile, exist=exist)
1025! IF (exist) THEN
1026! CALL LoadEsriAsciiGrid(Path, buildingsname, xllcorner, yllcorner, cellsize, NoData)
1027! ALLOCATE (buildings(1, 1))
1028! buildings = tempgrid
1029! DEALLOCATE (tempgrid)
1030! ELSE
1031! !!! Not ready, should return error. Also for the other grids
1032! ENDIF
1033
1034! ! Time related info
1035! !timestepdec=t_INTERVAL/(1440.*60.)
1036! timeadd = 0.00
1037
1038! ! Initiate map for surface temperature delay
1039! ALLOCATE (Tgmap1(1, 1))
1040! Tgmap1 = 0.0
1041
1042! RETURN
1043
1044! 274 CALL ErrorHint(40, TRIM("SOLWEIGinput.nml FileCode is missing"), notUsed, notUsed, notUsedI)
1045
1046! END SUBROUTINE SOLWEIG_Initial
1047
1048 SUBROUTINE kside_veg_v24( &
1049 shadow, F_sh, &
1050 radI, radG, radD, azimuth, altitude, psi, t, albedo, & ! input
1051 Keast, Knorth, Ksouth, Kwest) ! output
1052
1053 IMPLICIT NONE
1054
1055 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
1056 REAL(KIND(1D0)), INTENT(in) :: radI
1057 REAL(KIND(1D0)), INTENT(in) :: radG
1058 REAL(KIND(1D0)), INTENT(in) :: radD
1059 REAL(KIND(1D0)), INTENT(in) :: azimuth
1060 REAL(KIND(1D0)), INTENT(in) :: altitude
1061 REAL(KIND(1D0)), INTENT(in) :: psi
1062 REAL(KIND(1D0)), INTENT(in) :: t
1063 REAL(KIND(1D0)), INTENT(in) :: albedo
1064
1065 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: shadow, F_sh
1066 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: Keast, Knorth, Ksouth, Kwest
1067
1068 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW
1069 REAL(KIND(1D0)), DIMENSION(1, 1) :: viktveg, viktwall
1070
1071 ! assuming the following SVF to ONE
1072 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfE = 1
1073 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfS = 1
1074 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfW = 1
1075 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfN = 1
1076 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfEveg = 1
1077 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfSveg = 1
1078 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfWveg = 1
1079 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfNveg = 1
1080
1081 ! Internal grids
1082 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: svfviktbuveg
1083
1084 ALLOCATE (svfviktbuveg(1, 1))
1085 ! New reflection equation 2012-05-25
1086 vikttot = 4.4897
1087 azie = azimuth + t
1088 azis = azimuth - 90 + t
1089 aziw = azimuth - 180 + t
1090 azin = azimuth - 270 + t
1091 ! sunw=cos(altitude*(pi/180)); ! anngle to walls
1092 !! Kside with weights
1093 CALL kvikt_veg(svfe, svfeveg, vikttot, viktveg, viktwall)
1094 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
1095 IF (azimuth > (360 - t) .OR. azimuth <= (180 - t)) THEN
1096 keast = radi*shadow*cos(altitude*(pi/180))*sin(azie*(pi/180)) + &
1097 radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1098 ELSE
1099 keast = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1100 END IF
1101
1102 CALL kvikt_veg(svfs, svfsveg, vikttot, viktveg, viktwall)
1103 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
1104 IF (azimuth > (90 - t) .AND. azimuth <= (270 - t)) THEN
1105 ksouth = radi*shadow*cos(altitude*(pi/180))*sin(azis*(pi/180)) + &
1106 radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1107 ELSE
1108 ksouth = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1109 END IF
1110
1111 CALL kvikt_veg(svfw, svfwveg, vikttot, viktveg, viktwall)
1112 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
1113 IF (azimuth > (180 - t) .AND. azimuth <= (360 - t)) THEN
1114 kwest = radi*shadow*cos(altitude*(pi/180))*sin(aziw*(pi/180)) + &
1115 radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1116 ELSE
1117 kwest = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1118 END IF
1119
1120 CALL kvikt_veg(svfn, svfnveg, vikttot, viktveg, viktwall)
1121 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
1122 IF (azimuth <= (90 - t) .OR. azimuth > (270 - t)) THEN
1123 knorth = radi*shadow*cos(altitude*(pi/180))*sin(azin*(pi/180)) + &
1124 radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1125 ELSE
1126 knorth = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh) !*sin(altitude*(pi/180));
1127 END IF
1128
1129 DEALLOCATE (svfviktbuveg)
1130 END SUBROUTINE kside_veg_v24
1131
1132 SUBROUTINE kvikt_veg(isvf, isvfveg, vikttot, & !input
1133 viktveg, viktwall) ! output
1134
1135 IMPLICIT NONE
1136 REAL(KIND(1D0)), INTENT(in) :: vikttot
1137 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: isvf
1138 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: isvfveg
1139 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: viktveg, viktwall
1140 REAL(KIND(1D0)), DIMENSION(1, 1) :: svfvegbu
1141
1142 !! Least
1143 viktwall = (vikttot &
1144 - (63.227*isvf**6 - 161.51*isvf**5 &
1145 + 156.91*isvf**4 - 70.424*isvf**3 &
1146 + 16.773*isvf**2 - 0.4863*isvf))/vikttot
1147
1148 svfvegbu = (isvfveg + isvf - 1) ! Vegetation plus buildings
1149 viktveg = (vikttot &
1150 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1151 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1152 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1153 viktveg = viktveg - viktwall
1154 END SUBROUTINE kvikt_veg
1155
1156 SUBROUTINE lside_veg_v2( &
1157 Ldown2d, Lup2d, &
1158 altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, zen, ldown, svfalfa, &
1159 Least, Lnorth, Lsouth, Lwest)
1160 ! This m-file is the current one that estimates L from the four cardinal points 20100414
1161 IMPLICIT NONE
1162
1163 REAL(KIND(1D0)), INTENT(in) :: altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, ldown, zen
1164
1165 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: svfalfa
1166 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: Ldown2d, Lup2d
1167 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: Least, Lnorth, Lsouth, Lwest
1168
1169 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW, c
1170
1171 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: svfalfaE, svfalfaS, svfalfaW, svfalfaN
1172 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: alfaB, betaB, betasun
1173 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: Lground, Lrefl, Lsky, Lsky_allsky, Lveg, Lwallsh, Lwallsun
1174 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: viktonlywall, viktaveg, svfvegbu
1175 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: oneminussvfE, oneminussvfS, oneminussvfW, oneminussvfN
1176 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
1177 INTEGER, PARAMETER :: SOLWEIG_ldown = 0 !! force to 0, TS 13 Dec 2019
1178
1179 REAL(KIND(1D0)), DIMENSION(1, 1) :: viktveg, viktsky, viktrefl, viktwall
1180 REAL(KIND(1D0)), DIMENSION(1, 1) :: F_sh
1181
1182 ! assuming all SVF to ONE, TS 14 Dec 2019
1183 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfE = 1
1184 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfS = 1
1185 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfW = 1
1186 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfN = 1
1187 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfEveg = 1
1188 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfSveg = 1
1189 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfWveg = 1
1190 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfNveg = 1
1191 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfEaveg = 1
1192 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfSaveg = 1
1193 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfWaveg = 1
1194 REAL(KIND(1D0)), DIMENSION(1, 1), PARAMETER :: svfNaveg = 1
1195
1196 ALLOCATE (oneminussvfe(1, 1))
1197 ALLOCATE (oneminussvfs(1, 1))
1198 ALLOCATE (oneminussvfw(1, 1))
1199 ALLOCATE (oneminussvfn(1, 1))
1200 ALLOCATE (svfalfae(1, 1))
1201 ALLOCATE (svfalfas(1, 1))
1202 ALLOCATE (svfalfaw(1, 1))
1203 ALLOCATE (svfalfan(1, 1))
1204 ALLOCATE (alfab(1, 1))
1205 ALLOCATE (betab(1, 1))
1206 ALLOCATE (betasun(1, 1))
1207 ALLOCATE (lground(1, 1))
1208 ALLOCATE (lrefl(1, 1))
1209 ALLOCATE (lsky(1, 1))
1210 ALLOCATE (lsky_allsky(1, 1))
1211 ALLOCATE (lveg(1, 1))
1212 ALLOCATE (lwallsh(1, 1))
1213 ALLOCATE (lwallsun(1, 1))
1214 ALLOCATE (viktonlywall(1, 1))
1215 ALLOCATE (viktaveg(1, 1))
1216 ALLOCATE (svfvegbu(1, 1))
1217
1218 ! IF (ALLOCATED(viktwall)) DEALLOCATE (viktwall); ALLOCATE (viktwall(1, 1))
1219 ! IF (ALLOCATED(viktsky)) DEALLOCATE (viktsky); ALLOCATE (viktsky(1, 1))
1220 ! IF (ALLOCATED(viktveg)) DEALLOCATE (viktveg); ALLOCATE (viktveg(1, 1))
1221 ! IF (ALLOCATED(viktrefl)) DEALLOCATE (viktrefl); ALLOCATE (viktrefl(1, 1))
1222
1223 !Building height angle from svf
1224 oneminussvfe = 1.-svfe; WHERE (oneminussvfe <= 0) oneminussvfe = 0.000000001 ! avoiding log(0)
1225 oneminussvfs = 1.-svfs; WHERE (oneminussvfs <= 0) oneminussvfs = 0.000000001 ! avoiding log(0)
1226 oneminussvfw = 1.-svfw; WHERE (oneminussvfw <= 0) oneminussvfw = 0.000000001 ! avoiding log(0)
1227 oneminussvfn = 1.-svfn; WHERE (oneminussvfn <= 0) oneminussvfn = 0.000000001 ! avoiding log(0)
1228 svfalfae = asin(exp((log(oneminussvfe))/2))
1229 svfalfas = asin(exp((log(oneminussvfs))/2))
1230 svfalfaw = asin(exp((log(oneminussvfw))/2))
1231 svfalfan = asin(exp((log(oneminussvfn))/2))
1232
1233 vikttot = 4.4897
1234 aziw = azimuth + t
1235 azin = azimuth - 90 + t
1236 azie = azimuth - 180 + t
1237 azis = azimuth - 270 + t
1238
1239 CALL cylindric_wedge(zen, svfalfa, f_sh) !Fraction shadow on building walls based on sun altitude and svf
1240 ! F_sh(isnan(F_sh))=0.5;
1241 f_sh = 2.*f_sh - 1. !(cylindric_wedge scaled 0-1)
1242
1243 IF (solweig_ldown == 1) THEN
1244 c = 1 - ci
1245 lsky_allsky = emis_sky*sbc*((ta + 273.15)**4)*(1 - c) + c*sbc*((ta + 273.15)**4)
1246 ELSE
1247 lsky_allsky = ldown
1248 END IF
1249
1250 !! Least
1251 CALL lvikt_veg(svfe, svfeveg, svfeaveg, vikttot, &
1252 viktveg, viktsky, viktrefl, viktwall)
1253
1254 IF (altitude > 0) THEN ! daytime
1255 alfab = atan(svfalfae)
1256 betab = atan(tan((svfalfae)*f_sh))
1257 betasun = ((alfab - betab)/2) + betab
1258 IF (azimuth > (180 - t) .AND. azimuth <= (360 - t)) THEN
1259 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azie*(pi/180)))**4)* &
1260 viktwall*(1 - f_sh)*cos(betasun)*0.5
1261 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1262 ELSE
1263 lwallsun = 0
1264 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1265 END IF
1266 ELSE !nighttime
1267 lwallsun = 0
1268 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1269 END IF
1270 lsky = ((svfe + svfeveg - 1)*lsky_allsky)*viktsky*0.5
1271 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1272 lground = lup2d*0.5
1273 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1274 least = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1275
1276 !! Lsouth
1277 CALL lvikt_veg(svfs, svfsveg, svfsaveg, vikttot, &
1278 viktveg, viktsky, viktrefl, viktwall)
1279
1280 IF (altitude > 0) THEN ! daytime
1281 alfab = atan(svfalfas)
1282 betab = atan(tan((svfalfas)*f_sh))
1283 betasun = ((alfab - betab)/2) + betab
1284 IF (azimuth <= (90 - t) .OR. azimuth > (270 - t)) THEN
1285 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azis*(pi/180)))**4)* &
1286 viktwall*(1 - f_sh)*cos(betasun)*0.5
1287 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1288 ELSE
1289 lwallsun = 0
1290 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1291 END IF
1292 ELSE !nighttime
1293 lwallsun = 0
1294 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1295 END IF
1296 lsky = ((svfs + svfsveg - 1)*lsky_allsky)*viktsky*0.5
1297 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1298 lground = lup2d*0.5
1299 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1300 lsouth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1301
1302 !! Lwest
1303 CALL lvikt_veg(svfw, svfwveg, svfwaveg, vikttot, &
1304 viktveg, viktsky, viktrefl, viktwall)
1305
1306 IF (altitude > 0) THEN ! daytime
1307 alfab = atan(svfalfaw)
1308 betab = atan(tan((svfalfaw)*f_sh))
1309 betasun = ((alfab - betab)/2) + betab
1310 IF (azimuth > (360 - t) .OR. azimuth <= (180 - t)) THEN
1311 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(aziw*(pi/180)))**4)* &
1312 viktwall*(1 - f_sh)*cos(betasun)*0.5
1313 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1314 ELSE
1315 lwallsun = 0
1316 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1317 END IF
1318 ELSE !nighttime
1319 lwallsun = 0
1320 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1321 END IF
1322 lsky = ((svfw + svfwveg - 1)*lsky_allsky)*viktsky*0.5
1323 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1324 lground = lup2d*0.5
1325 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1326 lwest = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1327
1328 !! Lnorth
1329 CALL lvikt_veg(svfn, svfnveg, svfnaveg, vikttot, &
1330 viktveg, viktsky, viktrefl, viktwall)
1331
1332 IF (altitude > 0) THEN ! daytime
1333 alfab = atan(svfalfan)
1334 betab = atan(tan((svfalfan)*f_sh))
1335 betasun = ((alfab - betab)/2) + betab
1336 IF (azimuth > (90 - t) .AND. azimuth <= (270 - t)) THEN
1337 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azin*(pi/180)))**4)* &
1338 viktwall*(1 - f_sh)*cos(betasun)*0.5
1339 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1340 ELSE
1341 lwallsun = 0
1342 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1343 END IF
1344 ELSE !nighttime
1345 lwallsun = 0
1346 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1347 END IF
1348 lsky = ((svfn + svfnveg - 1)*lsky_allsky)*viktsky*0.5
1349 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1350 lground = lup2d*0.5
1351 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1352 lnorth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1353
1354 DEALLOCATE (svfalfae)
1355 DEALLOCATE (svfalfas)
1356 DEALLOCATE (svfalfaw)
1357 DEALLOCATE (svfalfan)
1358 DEALLOCATE (alfab)
1359 DEALLOCATE (betab)
1360 DEALLOCATE (betasun)
1361 DEALLOCATE (lground)
1362 DEALLOCATE (lrefl)
1363 DEALLOCATE (lsky)
1364 DEALLOCATE (lsky_allsky)
1365 DEALLOCATE (lveg)
1366 DEALLOCATE (lwallsh)
1367 DEALLOCATE (lwallsun)
1368 DEALLOCATE (viktonlywall)
1369 DEALLOCATE (viktaveg)
1370 DEALLOCATE (svfvegbu)
1371
1372 END SUBROUTINE lside_veg_v2
1373
1374 SUBROUTINE lvikt_veg(isvf, isvfveg, isvfaveg, vikttot, & ! input
1375 viktveg, viktsky, viktrefl, viktwall) !output
1376
1377 IMPLICIT NONE
1378 REAL(KIND(1D0)), INTENT(in) :: vikttot
1379 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: isvf
1380 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: isvfveg
1381 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: isvfaveg
1382 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: viktveg
1383 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: viktsky
1384 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: viktrefl
1385 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: viktwall
1386
1387 REAL(KIND(1D0)), DIMENSION(1, 1) :: viktonlywall
1388 REAL(KIND(1D0)), DIMENSION(1, 1) :: viktaveg
1389 REAL(KIND(1D0)), DIMENSION(1, 1) :: svfvegbu
1390
1391 !allocate(svfalfaE(1,1))
1392 !allocate(svfalfaS(1,1))
1393 !allocate(svfalfaW(1,1))
1394 !allocate(svfalfaN(1,1))
1395 !allocate(alfaB(1,1))
1396 !allocate(betaB(1,1))
1397
1398 !! Least
1399 viktonlywall = (vikttot - &
1400 (63.227*isvf**6 - 161.51*isvf**5 + 156.91*isvf**4 &
1401 - 70.424*isvf**3 + 16.773*isvf**2 - 0.4863*isvf))/vikttot
1402
1403 viktaveg = (vikttot &
1404 - (63.227*isvfaveg**6 - 161.51*isvfaveg**5 &
1405 + 156.91*isvfaveg**4 - 70.424*isvfaveg**3 &
1406 + 16.773*isvfaveg**2 - 0.4863*isvfaveg))/vikttot
1407
1408 viktwall = viktonlywall - viktaveg
1409
1410 svfvegbu = (isvfveg + isvf - 1) ! Vegetation plus buildings
1411 viktsky = (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1412 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1413 + 16.773*svfvegbu**2 - 0.4863*svfvegbu)/vikttot
1414 viktrefl = (vikttot &
1415 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1416 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1417 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1418 viktveg = (vikttot &
1419 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1420 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1421 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1422 viktveg = viktveg - viktwall
1423
1424 END SUBROUTINE lvikt_veg
1425 ! FUNCTION TO RETURN 0 IF IX=0, 1 IF 0<IX<MAXPOS+1,-1 OTHERWISE.
1426 ! MAXPOS is given as the maximum positive Integer.
1427 SUBROUTINE issign(IX, MAXPOS, ISIGNM)
1428 REAL(KIND(1D0)) IX, MAXPOS, ISIGNM
1429 isignm = 1.0
1430 IF (ix < 0 .OR. ix > maxpos) isignm = -1
1431 IF (ix == 0) isignm = 0
1432 RETURN
1433 END SUBROUTINE issign
1434
1435 !=====================================================
1436 ! fuction to convert interger to string
1437 CHARACTER(len=20) FUNCTION str(k)
1438 INTEGER, INTENT(in) :: k
1439 WRITE (str, *) k
1440 str = adjustl(str)
1441 END FUNCTION str
1442
1443 ! SUBROUTINE SaveGrids
1444
1445 ! IMPLICIT NONE
1446 ! CHARACTER(len=100) :: GridPath, GridName, GridPath2
1447 ! CHARACTER(len=4) ::doy, hour
1448 ! !real(kind(1d0)), allocatable, dimension(:,:):: savegrid
1449
1450 ! ALLOCATE (savegrid(1, 1))
1451
1452 ! IF (Tmrt_out == 1) THEN
1453 ! Gridpath2 = 'Grids/'
1454 ! GridPath = TRIM(FileOutputPath)//TRIM(GridPath2)
1455 ! WRITE (doy, '(i3)') id
1456 ! WRITE (hour, '(i2)') it
1457 ! hour = ADJUSTL(hour)
1458 ! doy = ADJUSTL(doy)
1459 ! GridName = 'Tmrt_'//TRIM(doy)//'_'//TRIM(hour)//'.txt'
1460 ! savegrid = Tmrt
1461 ! CALL SaveEsriAsciiGrid(GridPath, GridName, xllcorner, yllcorner, cellsize, NoData)
1462 ! END IF
1463 ! IF (Lup2d_out == 1) THEN
1464 ! Gridpath2 = 'Grids/'
1465 ! GridPath = TRIM(FileOutputPath)//TRIM(GridPath2)
1466 ! WRITE (doy, '(i3)') id
1467 ! WRITE (hour, '(i2)') it
1468 ! hour = ADJUSTL(hour)
1469 ! doy = ADJUSTL(doy)
1470 ! GridName = 'Tmrt_'//TRIM(doy)//'_'//TRIM(hour)//'.txt'
1471 ! savegrid = Tmrt
1472 ! CALL SaveEsriAsciiGrid(GridPath, GridName, xllcorner, yllcorner, cellsize, NoData)
1473 ! END IF
1474 ! IF (Ldown2d_out == 1) THEN
1475 ! Gridpath2 = 'Grids/'
1476 ! GridPath = TRIM(FileOutputPath)//TRIM(GridPath2)
1477 ! WRITE (doy, '(i3)') id
1478 ! WRITE (hour, '(i2)') it
1479 ! hour = ADJUSTL(hour)
1480 ! doy = ADJUSTL(doy)
1481 ! GridName = 'Ldown2d_'//TRIM(doy)//'_'//TRIM(hour)//'.txt'
1482 ! savegrid = Ldown2d
1483 ! CALL SaveEsriAsciiGrid(GridPath, GridName, xllcorner, yllcorner, cellsize, NoData)
1484 ! END IF
1485 ! IF (Kup2d_out == 1) THEN
1486 ! Gridpath2 = 'Grids/'
1487 ! GridPath = TRIM(FileOutputPath)//TRIM(GridPath2)
1488 ! WRITE (doy, '(i3)') id
1489 ! WRITE (hour, '(i2)') it
1490 ! hour = ADJUSTL(hour)
1491 ! doy = ADJUSTL(doy)
1492 ! GridName = 'Kup2d_'//TRIM(doy)//'_'//TRIM(hour)//'.txt'
1493 ! savegrid = Kup2d
1494 ! CALL SaveEsriAsciiGrid(GridPath, GridName, xllcorner, yllcorner, cellsize, NoData)
1495 ! END IF
1496 ! IF (Kdown2d_out == 1) THEN
1497 ! Gridpath2 = 'Grids/'
1498 ! GridPath = TRIM(FileOutputPath)//TRIM(GridPath2)
1499 ! WRITE (doy, '(i3)') id
1500 ! WRITE (hour, '(i2)') it
1501 ! hour = ADJUSTL(hour)
1502 ! doy = ADJUSTL(doy)
1503 ! GridName = 'Kdown2d_'//TRIM(doy)//'_'//TRIM(hour)//'.txt'
1504 ! savegrid = Kdown2d
1505 ! CALL SaveEsriAsciiGrid(GridPath, GridName, xllcorner, yllcorner, cellsize, NoData)
1506 ! END IF
1507 ! IF (GVF_out == 1) THEN
1508 ! Gridpath2 = 'Grids/'
1509 ! GridPath = TRIM(FileOutputPath)//TRIM(GridPath2)
1510 ! WRITE (doy, '(i3)') id
1511 ! WRITE (hour, '(i2)') it
1512 ! hour = ADJUSTL(hour)
1513 ! doy = ADJUSTL(doy)
1514 ! GridName = 'GVF_'//TRIM(doy)//'_'//TRIM(hour)//'.txt'
1515 ! savegrid = gvf
1516 ! CALL SaveEsriAsciiGrid(GridPath, GridName, xllcorner, yllcorner, cellsize, NoData)
1517 ! END IF
1518
1519 ! DEALLOCATE (savegrid)
1520
1521 ! END SUBROUTINE SaveGrids
1522 !------------------------------------------------!
1523 ! Shadow casting algorithm, ground and buildings !
1524 !------------------------------------------------!
1525 SUBROUTINE shadowingfunction_urban(azimuth, altitude, scale, shadow)
1526 ! This m.file calculates shadows on a DSM
1527 ! Code originate from Carlo Rattis thesis
1528 ! This code is translated from Matlab by
1529 ! Fredrik Lindberg, Gothenburg University
1530 ! Last modified LJ 27 Jan 2016 - Removal of tabs and fixing real-int conversions
1531
1532 IMPLICIT NONE
1533 REAL(KIND(1D0)), INTENT(in) :: azimuth, altitude, scale
1534 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: shadow
1535
1536 REAL(KIND(1D0)), DIMENSION(1, 1) :: a
1537
1538 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
1539 REAL(KIND(1D0)), PARAMETER :: maxpos = 10000000000.0
1540 !real, allocatable, dimension(:,:) :: a,f,temp !! Defined in matsize
1541 REAL(KIND(1D0)) :: degrees, azi, alt, dx, dy, dz, ds, absdx, absdy
1542 REAL(KIND(1D0)) :: amaxvalue, pibyfour, threetimespibyfour, fivetimespibyfour
1543 REAL(KIND(1D0)) :: seventimespibyfour, sinazimuth, cosazimuth, tanazimuth
1544 REAL(KIND(1D0)) :: signsinazimuth, signcosazimuth, dssin, dscos, tanaltitudebyscale
1545 INTEGER :: index, xc1, xc2, yc1, yc2, xp1, xp2, yp1, yp2 !,row,col !,1,1
1546 ! Internal grids
1547 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: temp, f
1548
1549 !special cases
1550 ! IF (altitude == 90) THEN
1551 ! altitude = altitude - 0.0001
1552 ! END IF
1553 ! IF (azimuth == 0) THEN
1554 ! azimuth = azimuth - 0.0001
1555 ! END IF
1556
1557 ! set a as ZERO for later initialisation
1558 a = 0
1559
1560 ! conversion
1561 degrees = pi/180
1562 azi = min(azimuth, 0 - 0.0001)*degrees
1563 alt = min(altitude, 90 - 0.0001)*degrees
1564
1565 ALLOCATE (f(1, 1))
1566 ALLOCATE (temp(1, 1))
1567
1568 ! IF (ALLOCATED(sh)) DEALLOCATE (sh)
1569 ! ALLOCATE (sh(1, 1))
1570
1571 ! initialise parameters
1572 f = a
1573 dx = 0
1574 dy = 0
1575 dz = 0
1576 temp = a*0.0
1577 index = 1
1578
1579 ! other loop parameters
1580 amaxvalue = maxval(a)
1581 pibyfour = pi/4.
1582 threetimespibyfour = 3.*pibyfour
1583 fivetimespibyfour = 5.*pibyfour
1584 seventimespibyfour = 7.*pibyfour
1585 sinazimuth = sin(azi)
1586 cosazimuth = cos(azi)
1587 tanazimuth = tan(azi)
1588 CALL issign(sinazimuth, maxpos, signsinazimuth)
1589 CALL issign(cosazimuth, maxpos, signcosazimuth)
1590 !signsinazimuth=sinazimuth/abs(sinazimuth)
1591 !signcosazimuth=cosazimuth/abs(cosazimuth)
1592 dssin = abs(1./sinazimuth)
1593 dscos = abs(1./cosazimuth)
1594 tanaltitudebyscale = tan(alt)/scale
1595
1596 ! main loop
1597 DO WHILE (amaxvalue >= dz .AND. abs(dx) <= 1 .AND. abs(dy) <= 1)
1598
1599 IF ((pibyfour <= azi .AND. azi < threetimespibyfour) .OR. (fivetimespibyfour <= azi .AND. azi < seventimespibyfour)) THEN
1600 dy = signsinazimuth*index
1601 dx = -1.*signcosazimuth*abs(nint(index/tanazimuth))
1602 ds = dssin
1603 ELSE
1604 dy = signsinazimuth*abs(nint(index*tanazimuth))
1605 dx = -1.*signcosazimuth*index
1606 ds = dscos
1607 END IF
1608
1609 dz = ds*index*tanaltitudebyscale
1610 temp = temp*0
1611
1612 absdx = abs(dx)
1613 absdy = abs(dy)
1614
1615 xc1 = int((dx + absdx)/2) + 1
1616 xc2 = (1 + int((dx - absdx)/2))
1617 yc1 = int((dy + absdy)/2) + 1
1618 yc2 = (1 + int((dy - absdy)/2))
1619 xp1 = -int((dx - absdx)/2) + 1
1620 xp2 = (1 - int((dx + absdx)/2))
1621 yp1 = -int((dy - absdy)/2) + 1
1622 yp2 = (1 - int((dy + absdy)/2))
1623
1624 temp(xp1:xp2, yp1:yp2) = a(xc1:xc2, yc1:yc2) - dz
1625
1626 f = max(f, temp)
1627 index = index + 1
1628
1629 END DO
1630
1631 f = f - a
1632 WHERE (f > 0)
1633 f = -1
1634 END WHERE
1635 shadow = f + 1
1636 !sh=f ! invert as in shadowingfunctionglobalradiation
1637 DEALLOCATE (f)
1638 DEALLOCATE (temp)
1639
1640 END SUBROUTINE shadowingfunction_urban
1641 !------------------------------------------------!
1642 ! Shadow casting algorithm, vegetation !
1643 !------------------------------------------------!
1644
1645 ! SUBROUTINE shadowingfunction_veg(azimuth, altitude, scale, amaxvalue)
1646
1647 ! ! This m.file calculates shadows on a DSM and for vegetation units
1648 ! ! This code is translated from Matlab by Fredrik Lindberg, Gothenburg University
1649
1650 ! IMPLICIT NONE
1651 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
1652 ! REAL(KIND(1d0)), PARAMETER :: maxpos = 10000000000.0
1653 ! REAL(KIND(1d0)) :: degrees, azi, alt, dx, dy, dz, ds, absdx, absdy, azimuth, altitude
1654 ! REAL(KIND(1d0)) :: amaxvalue, pibyfour, threetimespibyfour, fivetimespibyfour
1655 ! REAL(KIND(1d0)) :: seventimespibyfour, sinazimuth, cosazimuth, tanazimuth
1656 ! REAL(KIND(1d0)) :: signsinazimuth, signcosazimuth, dssin, dscos, tanaltitudebyscale, scale
1657 ! INTEGER :: index, xc1, xc2, yc1, yc2, xp1, xp2, yp1, yp2!,test!,row,col !,1,1
1658 ! ! Internal grids
1659 ! REAL(KIND(1d0)), ALLOCATABLE, DIMENSION(:, :) :: f, temp, tmp, stopbuild, stopveg, g, bushplant, tempvegdem, tempvegdem2
1660 ! REAL(KIND(1d0)), ALLOCATABLE, DIMENSION(:, :) :: fabovea, gabovea, tempbush, firstvegdem, vegsh2
1661 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: bush, vegdem, vegdem2
1662 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: a, sh, vbshvegsh, vegsh
1663
1664 ! !real :: start_time,end_time
1665
1666 ! !special case
1667 ! IF (altitude == 90) THEN
1668 ! altitude = altitude - 0.0001
1669 ! END IF
1670 ! IF (azimuth == 0) THEN
1671 ! azimuth = azimuth - 0.0001
1672 ! END IF
1673
1674 ! ! conversion
1675 ! degrees = pi/180;
1676 ! azi = azimuth*degrees;
1677 ! alt = altitude*degrees;
1678 ! ! IF (ALLOCATED(sh)) DEALLOCATE (sh)
1679 ! ! ALLOCATE (sh(1, 1))
1680 ! ! IF (ALLOCATED(vegsh)) DEALLOCATE (vegsh)
1681 ! ! ALLOCATE (vegsh(1, 1))
1682 ! ! IF (ALLOCATED(vbshvegsh)) DEALLOCATE (vbshvegsh)
1683 ! ! ALLOCATE (vbshvegsh(1, 1))
1684
1685 ! ! allocation of grids
1686 ! ALLOCATE (f(1, 1))
1687 ! ALLOCATE (temp(1, 1))
1688 ! ALLOCATE (tmp(1, 1))
1689 ! ALLOCATE (stopbuild(1, 1))
1690 ! ALLOCATE (stopveg(1, 1))
1691 ! ALLOCATE (g(1, 1))
1692 ! ALLOCATE (bushplant(1, 1))
1693 ! ALLOCATE (tempvegdem(1, 1))
1694 ! ALLOCATE (tempvegdem2(1, 1))
1695 ! ALLOCATE (fabovea(1, 1))
1696 ! ALLOCATE (gabovea(1, 1))
1697 ! ALLOCATE (firstvegdem(1, 1))
1698 ! ALLOCATE (tempbush(1, 1))
1699 ! ALLOCATE (vegsh2(1, 1))
1700
1701 ! ! initialise parameters
1702 ! f = a
1703 ! dx = 0
1704 ! dy = 0
1705 ! dz = 0
1706 ! temp = a*0.0
1707 ! sh = temp
1708 ! vegsh = sh
1709 ! stopbuild = sh
1710 ! stopveg = sh
1711 ! vbshvegsh = sh
1712 ! g = sh
1713 ! bushplant = temp
1714 ! WHERE (bush > 1)
1715 ! bushplant = 1
1716 ! END WHERE
1717
1718 ! index = 1
1719 ! !test=0
1720 ! ! other loop parameters
1721 ! !amaxvalue=maxval(a)
1722 ! pibyfour = pi/4.
1723 ! threetimespibyfour = 3.*pibyfour;
1724 ! fivetimespibyfour = 5.*pibyfour;
1725 ! seventimespibyfour = 7.*pibyfour;
1726 ! sinazimuth = SIN(azi);
1727 ! cosazimuth = COS(azi);
1728 ! tanazimuth = TAN(azi);
1729 ! CALL issign(sinazimuth, maxpos, signsinazimuth)
1730 ! CALL issign(cosazimuth, maxpos, signcosazimuth)
1731 ! !signsinazimuth=sinazimuth/abs(sinazimuth);
1732 ! !signcosazimuth=cosazimuth/abs(cosazimuth);
1733 ! dssin = ABS(1./sinazimuth);
1734 ! dscos = ABS(1./cosazimuth);
1735 ! tanaltitudebyscale = TAN(alt)/scale;
1736 ! DO WHILE (amaxvalue >= dz .AND. ABS(dx) <= 1 .AND. ABS(dy) <= 1)
1737
1738 ! IF ((pibyfour <= azi .AND. azi < threetimespibyfour) .OR. (fivetimespibyfour <= azi .AND. azi < seventimespibyfour)) THEN
1739 ! dy = signsinazimuth*index
1740 ! dx = -1.*signcosazimuth*ABS(NINT(index/tanazimuth))
1741 ! ds = dssin
1742 ! ELSE
1743 ! dy = signsinazimuth*ABS(NINT(index*tanazimuth))
1744 ! dx = -1.*signcosazimuth*index
1745 ! ds = dscos
1746 ! END IF
1747
1748 ! dz = ds*index*tanaltitudebyscale
1749 ! temp = temp*0
1750 ! tempvegdem = temp
1751 ! tempvegdem2 = temp
1752
1753 ! absdx = ABS(dx)
1754 ! absdy = ABS(dy)
1755
1756 ! xc1 = INT(((dx + absdx)/2)) + 1
1757 ! xc2 = (1 + INT((dx - absdx)/2))
1758 ! yc1 = INT((dy + absdy)/2) + 1
1759 ! yc2 = (1 + INT((dy - absdy)/2))
1760 ! xp1 = -INT((dx - absdx)/2) + 1
1761 ! xp2 = (1 - INT((dx + absdx)/2))
1762 ! yp1 = -INT((dy - absdy)/2) + 1
1763 ! yp2 = (1 - INT((dy + absdy)/2))
1764
1765 ! temp(xp1:xp2, yp1:yp2) = a(xc1:xc2, yc1:yc2) - dz
1766 ! tempvegdem(xp1:xp2, yp1:yp2) = vegdem(xc1:xc2, yc1:yc2) - dz
1767 ! tempvegdem2(xp1:xp2, yp1:yp2) = vegdem2(xc1:xc2, yc1:yc2) - dz
1768
1769 ! f = MAX(f, temp)
1770 ! WHERE (f > a) !sh(f>a)=1;sh(f<=a)=0; !Moving building shadow
1771 ! sh = 1
1772 ! ELSEWHERE
1773 ! sh = 0
1774 ! END WHERE
1775 ! WHERE (tempvegdem > a) !fabovea=tempvegdem>a; !vegdem above DEM
1776 ! fabovea = 1
1777 ! ELSEWHERE
1778 ! fabovea = 0
1779 ! END WHERE
1780 ! WHERE (tempvegdem2 > a) !gabovea=tempvegdem2>a; !vegdem2 above DEM
1781 ! gabovea = 1
1782 ! ELSEWHERE
1783 ! gabovea = 0
1784 ! END WHERE
1785 ! vegsh2 = fabovea - gabovea
1786 ! vegsh = MAX(vegsh, vegsh2)
1787 ! WHERE ((vegsh*sh) > 0) !vegsh(vegsh.*sh>0)=0;! removing shadows 'behind' buildings
1788 ! vegsh = 0
1789 ! END WHERE
1790 ! vbshvegsh = vegsh + vbshvegsh
1791
1792 ! ! vegsh at high sun altitudes
1793 ! IF (index == 1) THEN
1794 ! firstvegdem = tempvegdem - temp
1795 ! WHERE (firstvegdem <= 0)!firstvegdem(firstvegdem<=0)=1000;
1796 ! firstvegdem = 1000
1797 ! END WHERE
1798 ! WHERE (firstvegdem < dz)!vegsh(firstvegdem<dz)=1;
1799 ! vegsh = 1
1800 ! END WHERE
1801 ! tmp = temp*0.0
1802 ! WHERE (vegdem2 > a)!vegsh=vegsh.*(vegdem2>a);
1803 ! tmp = 1
1804 ! END WHERE
1805 ! vegsh = vegsh*tmp
1806 ! vbshvegsh = temp*0.0 !vbshvegsh=zeros(1,1);
1807 ! END IF
1808
1809 ! ! Bush shadow on bush plant
1810 ! tmp = fabovea*bush
1811 ! IF ((MAXVAL(bush) > 0) .AND. (MAXVAL(tmp) > 0)) THEN
1812 ! tempbush = temp*0.0
1813 ! tempbush(xp1:xp2, yp1:yp2) = bush(xc1:xc2, yc1:yc2) - dz
1814 ! g = MAX(g, tempbush)
1815 ! g = bushplant*g
1816 ! END IF
1817
1818 ! index = index + 1
1819 ! END DO
1820
1821 ! sh = 1 - sh
1822 ! WHERE (vbshvegsh > 0)!vbshvegsh(vbshvegsh>0)=1;
1823 ! vbshvegsh = 1
1824 ! END WHERE
1825 ! vbshvegsh = vbshvegsh - vegsh;
1826 ! IF (MAXVAL(bush) > 0) THEN
1827 ! g = g - bush
1828 ! WHERE (g > 0)!g(g>0)=1;g(g<0)=0;
1829 ! g = 1
1830 ! ELSEWHERE
1831 ! g = 0
1832 ! END WHERE
1833 ! vegsh = vegsh - bushplant + g
1834 ! WHERE (vegsh < 0)!vegsh(vegsh<0)=0;
1835 ! vegsh = 0
1836 ! END WHERE
1837 ! END IF
1838
1839 ! WHERE (vegsh > 0)!vegsh(vegsh>0)=1;
1840 ! vegsh = 1
1841 ! END WHERE
1842 ! vegsh = 1 - vegsh
1843 ! vbshvegsh = 1 - vbshvegsh
1844
1845 ! !deallocation of grids
1846 ! DEALLOCATE (f)
1847 ! DEALLOCATE (temp)
1848 ! DEALLOCATE (tmp)
1849 ! DEALLOCATE (stopbuild)
1850 ! DEALLOCATE (stopveg)
1851 ! DEALLOCATE (g)
1852 ! DEALLOCATE (bushplant)
1853 ! DEALLOCATE (tempvegdem)
1854 ! DEALLOCATE (tempvegdem2)
1855 ! DEALLOCATE (fabovea)
1856 ! DEALLOCATE (gabovea)
1857 ! DEALLOCATE (firstvegdem)
1858 ! DEALLOCATE (tempbush)
1859 ! DEALLOCATE (vegsh2)
1860
1861 ! END SUBROUTINE shadowingfunction_veg
1862
1863 SUBROUTINE sunonsurface_veg(iazimuthA, scale, buildings, first, second, psi, sos)
1864 ! This m-file creates a boolean image of sunlit walls.
1865 ! Shadows from both buildings and vegetation is accounted for
1866 ! moving building in the direction of the sun
1867 ! Last modified:
1868 ! LJ 27 Jan 2016 - Removal of tabs and fixing real-int conversions
1869 IMPLICIT NONE
1870
1871 INTEGER :: first, second
1872
1873 REAL(KIND(1D0)), INTENT(in) :: iazimuthA
1874 REAL(KIND(1D0)), INTENT(in) :: scale
1875 REAL(KIND(1D0)), INTENT(in) :: psi
1876
1877 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(in) :: buildings
1878 REAL(KIND(1D0)), DIMENSION(1, 1), INTENT(out) :: sos
1879
1880 REAL(KIND(1D0)), DIMENSION(1, 1) :: sunwall
1881
1882 REAL(KIND(1D0)), DIMENSION(1, 1) :: sh, vegsh
1883
1884 REAL(KIND(1D0)) :: iazimuth, sinazimuth, cosazimuth, tanazimuth
1885 INTEGER :: index, xc1, xc2, yc1, yc2, xp1, xp2, yp1, yp2, n
1886 REAL(KIND(1D0)) :: dx, dy, ds, absdx, absdy !,dz
1887 REAL(KIND(1D0)) :: pibyfour, threetimespibyfour, fivetimespibyfour
1888 REAL(KIND(1D0)) :: seventimespibyfour
1889 REAL(KIND(1D0)) :: signsinazimuth, signcosazimuth, dssin, dscos
1890 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: weightsumwall, weightsumsh, gvf1, gvf2
1891 REAL(KIND(1D0)), ALLOCATABLE, DIMENSION(:, :) :: f, tempsh, tempbu, tempbub, tempwallsun, tempb, sh1
1892
1893 REAL(KIND(1D0)), PARAMETER :: pi = 3.141592653589793
1894 REAL(KIND(1D0)), PARAMETER :: maxpos = 10000000000.0
1895
1896 ALLOCATE (weightsumwall(1, 1))
1897 ALLOCATE (weightsumsh(1, 1))
1898 ALLOCATE (gvf1(1, 1))
1899 ALLOCATE (gvf2(1, 1))
1900 ALLOCATE (f(1, 1))
1901 ALLOCATE (tempsh(1, 1))
1902 ALLOCATE (tempbu(1, 1))
1903 ALLOCATE (tempbub(1, 1))
1904 ALLOCATE (tempwallsun(1, 1))
1905 ALLOCATE (tempb(1, 1))
1906 ALLOCATE (sh1(1, 1))
1907
1908 iazimuth = iazimutha*(pi/180)
1909 !special cases
1910 IF (iazimuth == 0) THEN
1911 iazimuth = iazimuth + 0.000001
1912 END IF
1913 ! loop parameters
1914 index = 0
1915 f = buildings
1916
1917 ! TODO: what are these:
1918 sh = 1
1919 vegsh = 1
1920 sunwall = 0
1921
1922 sh1 = sh - (1 - vegsh)*(1 - psi)
1923 dx = 0
1924 dy = 0
1925 ds = 0
1926
1927 tempsh = 0.0d0
1928 tempbu = 0.0d0
1929 tempbub = 0.0d0
1930 tempwallsun = 0.0d0
1931 !sh = 0.0D0
1932 weightsumsh = 0.0d0
1933 weightsumwall = 0.0d0
1934
1935 first = int(real(first, kind(1d0))*scale) !Int added around the equation as first and second are both integers
1936 second = int(real(second, kind(1d0))*scale)
1937
1938 ! other loop parameters
1939 pibyfour = pi/4.
1940 threetimespibyfour = 3.*pibyfour
1941 fivetimespibyfour = 5.*pibyfour
1942 seventimespibyfour = 7.*pibyfour
1943 sinazimuth = sin(iazimuth)
1944 cosazimuth = cos(iazimuth)
1945 tanazimuth = tan(iazimuth)
1946 CALL issign(sinazimuth, maxpos, signsinazimuth)
1947 CALL issign(cosazimuth, maxpos, signcosazimuth)
1948 !signsinazimuth=sinazimuth/abs(sinazimuth)
1949 !signcosazimuth=cosazimuth/abs(cosazimuth)
1950 dssin = abs(1./sinazimuth)
1951 dscos = abs(1./cosazimuth)
1952
1953 !! The Shadow casting algorithm
1954 DO n = 1, second
1955 IF ((pibyfour <= iazimuth .AND. iazimuth < threetimespibyfour) &
1956 .OR. (fivetimespibyfour <= iazimuth .AND. iazimuth < seventimespibyfour)) THEN
1957 dy = signsinazimuth*index
1958 dx = -1.*signcosazimuth*abs(nint(index/tanazimuth))
1959 ds = dssin
1960 ELSE
1961 dy = signsinazimuth*abs(nint(index*tanazimuth))
1962 dx = -1.*signcosazimuth*index
1963 ds = dscos
1964 END IF
1965
1966 absdx = abs(dx)
1967 absdy = abs(dy)
1968
1969 xc1 = int((dx + absdx)/2) + 1
1970 xc2 = (1 + int((dx - absdx)/2))
1971 yc1 = int((dy + absdy)/2) + 1
1972 yc2 = (1 + int((dy - absdy)/2))
1973 xp1 = -int((dx - absdx)/2) + 1
1974 xp2 = (1 - int((dx + absdx)/2))
1975 yp1 = -int((dy - absdy)/2) + 1
1976 yp2 = (1 - int((dy + absdy)/2))
1977
1978 tempbu(xp1:xp2, yp1:yp2) = buildings(xc1:xc2, yc1:yc2) !moving building
1979
1980 tempsh(xp1:xp2, yp1:yp2) = sh1(xc1:xc2, yc1:yc2) !moving shadow image
1981 f = min(f, tempbu) !utsmetning of buildings
1982
1983 weightsumsh = weightsumsh + tempsh*f
1984
1985 tempwallsun(xp1:xp2, yp1:yp2) = sunwall(xc1:xc2, yc1:yc2) !moving building wall in sun image
1986 tempb = tempwallsun*f
1987 WHERE ((tempb + tempbub) > 0) !tempbub=(tempb+tempbub)>0==1
1988 tempbub = 1.
1989 END WHERE
1990
1991 weightsumwall = weightsumwall + tempbub
1992
1993 IF (index*scale == first) THEN
1994 gvf1 = (weightsumwall + weightsumsh)/first
1995 WHERE (gvf1 > 1)
1996 gvf1 = 1.
1997 END WHERE
1998 END IF
1999 index = index + 1
2000
2001 END DO
2002 gvf2 = (weightsumsh + weightsumwall)/second
2003 WHERE (gvf2 > 1)
2004 gvf2 = 1.
2005 END WHERE
2006
2007 ! Weighting
2008 sos = (gvf1*0.5 + gvf2*0.4)/0.9
2009
2010 DEALLOCATE (weightsumwall)
2011 DEALLOCATE (weightsumsh)
2012 DEALLOCATE (gvf1)
2013 DEALLOCATE (gvf2)
2014 DEALLOCATE (f)
2015 DEALLOCATE (tempsh)
2016 DEALLOCATE (tempbu)
2017 DEALLOCATE (tempbub)
2018 DEALLOCATE (tempwallsun)
2019 DEALLOCATE (tempb)
2020 DEALLOCATE (sh1)
2021
2022 END SUBROUTINE sunonsurface_veg
2023
2025 ! SUBROUTINE wallinsun_veg(azimuth,sunwall)
2026 ! ! This m-file creates a boolean image of sunlit walls.
2027 ! ! Shadows from both buildings and vegetation is accounted for
2028 ! ! moving building in the direction of the sun
2029 ! ! Last modified:
2030 ! ! LJ 27 Jan 2017 - Change of equations xc1...yp2 to account for the change from real to integer
2031 ! !---------------------------------------------------------------------------------
2032
2033 ! IMPLICIT NONE
2034 ! REAL(KIND(1d0)) ,intent(in) :: azimuth
2035 ! REAL(KIND(1d0)), DIMENSION(1, 1),intent(out) :: sunwall
2036
2037 ! REAL(KIND(1d0)) ::iazimuth
2038 ! INTEGER :: index, xc1, xc2, yc1, yc2, xp1, xp2, yp1, yp2
2039 ! REAL(KIND(1d0)) :: dx, dy, dz, ds, absdx, absdy
2040 ! REAL(KIND(1d0)) :: pibyfour, threetimespibyfour, fivetimespibyfour
2041 ! REAL(KIND(1d0)) :: seventimespibyfour, sinazimuth, cosazimuth, tanazimuth
2042 ! REAL(KIND(1d0)) :: signsinazimuth, signcosazimuth, dssin, dscos, azi
2043 ! REAL(KIND(1d0)), PARAMETER :: pi = 3.141592653589793
2044 ! REAL(KIND(1d0)), PARAMETER :: maxpos = 10000000000.0
2045 ! ! Internal grids
2046 ! REAL(KIND(1d0)), ALLOCATABLE, DIMENSION(:, :) :: temp, sh1
2047 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: sh, vegsh
2048 ! REAL(KIND(1d0)), DIMENSION(1, 1) :: buildings
2049
2050 ! ALLOCATE (temp(1, 1))
2051 ! ALLOCATE (sh1(1, 1))
2052 ! !allocate(vegsh(1,1))
2053 ! !allocate(a(1,1))
2054 ! !allocate(buildings(1,1))
2055
2056 ! ! IF (ALLOCATED(sunwall)) DEALLOCATE (sunwall); ALLOCATE (sunwall(1, 1))
2057
2058 ! iazimuth = azimuth + 180
2059 ! IF (iazimuth >= 360) THEN
2060 ! iazimuth = iazimuth - 360
2061 ! END IF
2062 ! !special cases
2063 ! IF (iazimuth == 0) THEN
2064 ! iazimuth = iazimuth + 0.00001
2065 ! END IF
2066 ! ! conversion into radians
2067 ! azi = iazimuth*(pi/180)
2068
2069 ! index = 1
2070 ! dx = 0.
2071 ! dy = 0.
2072 ! dz = 0.
2073 ! ds = 0.
2074
2075 ! temp = 0.0D0
2076
2077 ! ! other loop parameters
2078 ! pibyfour = pi/4.
2079 ! threetimespibyfour = 3.*pibyfour
2080 ! fivetimespibyfour = 5.*pibyfour
2081 ! seventimespibyfour = 7.*pibyfour
2082 ! sinazimuth = SIN(azi)
2083 ! cosazimuth = COS(azi)
2084 ! tanazimuth = TAN(azi)
2085 ! CALL issign(sinazimuth, maxpos, signsinazimuth)
2086 ! CALL issign(cosazimuth, maxpos, signcosazimuth)
2087 ! !signsinazimuth=sinazimuth/abs(sinazimuth)
2088 ! !signcosazimuth=cosazimuth/abs(cosazimuth)
2089 ! dssin = ABS(1./sinazimuth)
2090 ! dscos = ABS(1./cosazimuth)
2091
2092 ! sh1 = vegsh + sh - 1.
2093 ! !! The Shadow casting algoritm
2094 ! IF ((pibyfour <= azi .AND. azi < threetimespibyfour) .OR. (fivetimespibyfour <= azi .AND. azi < seventimespibyfour)) THEN
2095 ! dy = signsinazimuth*index
2096 ! dx = -1.*signcosazimuth*ABS(NINT(index/tanazimuth))
2097 ! ds = dssin
2098 ! ELSE
2099 ! dy = signsinazimuth*ABS(NINT(index*tanazimuth))
2100 ! dx = -1.*signcosazimuth*index
2101 ! ds = dscos
2102 ! END IF
2103
2104 ! ! note: dx and dy represent absolute values while ds is an incremental value
2105
2106 ! absdx = ABS(dx)
2107 ! absdy = ABS(dy)
2108
2109 ! xc1 = INT((dx + absdx)/2) + 1 !LJ added int to the equation to account for the conversion from real to int
2110 ! xc2 = (1 + INT((dx - absdx)/2))
2111 ! yc1 = INT((dy + absdy)/2) + 1
2112 ! yc2 = (1 + INT((dy - absdy)/2))
2113 ! xp1 = -INT((dx - absdx)/2) + 1
2114 ! xp2 = (1 - INT((dx + absdx)/2))
2115 ! yp1 = -INT((dy - absdy)/2) + 1
2116 ! yp2 = (1 - INT((dy + absdy)/2))
2117
2118 ! temp(xp1:xp2, yp1:yp2) = buildings(xc1:xc2, yc1:yc2)
2119
2120 ! sunwall = temp - buildings
2121 ! WHERE (sunwall == 1) !f1(f1==1)=0
2122 ! sunwall = 0
2123 ! END WHERE
2124 ! WHERE (sunwall == -1) !f1(f1==-1)=1
2125 ! sunwall = 1
2126 ! END WHERE
2127 ! sunwall = sh1*sunwall
2128
2129 ! DEALLOCATE (temp)
2130 ! DEALLOCATE (sh1)
2131
2132 ! END SUBROUTINE wallinsun_veg
2133
2134END MODULE solweig_module
allocatearray::a2
real(kind(1d0)) a2
Definition suews_ctrl_const.f95:420
allocatearray::deg2rad
real(kind(1d0)), parameter deg2rad
Definition suews_ctrl_const.f95:412
allocatearray::ncolumnsdataoutsolweig
integer, parameter ncolumnsdataoutsolweig
Definition suews_ctrl_const.f95:80
allocatearray::rad2deg
real(kind(1d0)), parameter rad2deg
Definition suews_ctrl_const.f95:412
defaultnotused::notused
real(kind(1d0)) notused
Definition suews_ctrl_const.f95:1394
solweig_module
Definition suews_phys_solweig.f95:1
solweig_module::hwtosvf_ground
real(kind(1d0)) function hwtosvf_ground(hw)
Definition suews_phys_solweig.f95:481
solweig_module::issign
subroutine issign(ix, maxpos, isignm)
Definition suews_phys_solweig.f95:1428
solweig_module::hwtosvf_roof
real(kind(1d0)) function hwtosvf_roof(hw)
Definition suews_phys_solweig.f95:496
solweig_module::cylindric_wedge
subroutine cylindric_wedge(zen, svfalfa, f_sh)
Definition suews_phys_solweig.f95:641
solweig_module::solweig_cal_main
subroutine solweig_cal_main(id, it, dectime, lamdap, lamdaf, avkdn, ldown, temp_c, avrh, press_hpa, tg, lat, zenith_deg, azimuth, scale, alb_ground, alb_bldg, emis_ground, emis_wall, heightgravity, dataoutlinesolweig)
Definition suews_phys_solweig.f95:75
solweig_module::clearnessindex_2013b
subroutine clearnessindex_2013b(zen, doy, ta, rh, radg, lat, p_kpa, i0, ci, kt, i0et, ciuncorr)
Definition suews_phys_solweig.f95:511
solweig_module::lvikt_veg
subroutine lvikt_veg(isvf, isvfveg, isvfaveg, vikttot, viktveg, viktsky, viktrefl, viktwall)
Definition suews_phys_solweig.f95:1376
solweig_module::sunonsurface_veg
subroutine sunonsurface_veg(iazimutha, scale, buildings, first, second, psi, sos)
Definition suews_phys_solweig.f95:1864
solweig_module::diffusefraction
subroutine diffusefraction(radg, altitude, kt, ta, rh, radi, radd)
Definition suews_phys_solweig.f95:724
solweig_module::shadowingfunction_urban
subroutine shadowingfunction_urban(azimuth, altitude, scale, shadow)
Definition suews_phys_solweig.f95:1526
solweig_module::kvikt_veg
subroutine kvikt_veg(isvf, isvfveg, vikttot, viktveg, viktwall)
Definition suews_phys_solweig.f95:1134
solweig_module::str
character(len=20) function str(k)
Definition suews_phys_solweig.f95:1438
solweig_module::kside_veg_v24
subroutine kside_veg_v24(shadow, f_sh, radi, radg, radd, azimuth, altitude, psi, t, albedo, keast, knorth, ksouth, kwest)
Definition suews_phys_solweig.f95:1052
solweig_module::cal_ratio_height2width
real(kind(1d0)) function cal_ratio_height2width(lamdap, lamdaf)
Definition suews_phys_solweig.f95:463
solweig_module::sun_distance
subroutine sun_distance(jday, d)
Definition suews_phys_solweig.f95:627
solweig_module::lside_veg_v2
subroutine lside_veg_v2(ldown2d, lup2d, altitude, ta, tw, sbc, emis_wall, emis_sky, t, ci, azimuth, zen, ldown, svfalfa, least, lnorth, lsouth, lwest)
Definition suews_phys_solweig.f95:1160
daylen
subroutine daylen(doy, xlat, dayl, dec, sndn, snup)
Definition suews_util_time.f95:227
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0