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