81 integer,
intent(in) ::id
82 integer,
intent(in) ::it
83 REAL(KIND(1d0)),
INTENT(in)::lamdaP
84 REAL(KIND(1d0)),
INTENT(in)::lamdaF
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
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
102 REAL(KIND(1D0)),
parameter:: absL = 0.97
103 REAL(KIND(1D0)),
parameter:: absK = 0.7
104 REAL(KIND(1D0)),
intent(in)::heightgravity
106 REAL(KIND(1D0)),
intent(in)::alb_bldg
107 REAL(KIND(1D0)),
intent(in)::alb_ground
108 REAL(KIND(1D0)),
intent(in)::emis_wall
109 REAL(KIND(1D0)),
intent(in)::emis_ground
111 REAL(KIND(1D0)),
DIMENSION(ncolumnsDataOutSol - 5),
INTENT(OUT) ::dataOutLineSOLWEIG
113 REAL(KIND(1d0)) :: t, Tstart, height, psi
114 REAL(KIND(1d0)) :: altitude, zen
115 REAL(KIND(1d0)) :: CI, c, I0, Kt, Tw, weight1
116 REAL(KIND(1d0)) :: Ta, RH, P, radG, radD, radI
117 REAL(KIND(1d0)) :: I0et, CIuncorr
118 REAL(KIND(1d0)) :: SNDN, SNUP, DEC, DAYL
119 REAL(KIND(1d0)) :: msteg, emis_sky, ea
121 INTEGER :: DOY, hour, first, second, j
122 REAL(KIND(1d0)) :: timeadd
123 REAL(KIND(1d0)) :: firstdaytime
125 REAL(KIND(1d0)) :: CIlatenight
126 REAL(KIND(1d0)) :: Fside
127 REAL(KIND(1d0)) :: Fup
128 REAL(KIND(1d0)) :: HW
130 REAL(KIND(1d0)),
DIMENSION(1, 1) :: svfalfa, sos, Tgmap1
131 REAL(KIND(1d0)),
DIMENSION(1, 1) :: 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
141 REAL(KIND(1d0)),
ALLOCATABLE,
DIMENSION(:, :) :: tmp, Knight, svf_blgd_veg, Tgmap0
143 REAL(KIND(1d0)),
PARAMETER :: azimuthA(1:18) = [(j*(360.0/18.0), j=0, 17)]
145 REAL(KIND(1d0)),
PARAMETER :: pi = 3.141592653589793
146 REAL(KIND(1d0)),
PARAMETER :: SBC = 5.67051e-8
149 INTEGER,
PARAMETER:: onlyglobal = 1
150 INTEGER,
PARAMETER:: usevegdem = 0
151 INTEGER,
PARAMETER:: row = 1
152 INTEGER,
PARAMETER:: col = 1
153 INTEGER,
PARAMETER:: Posture = 1
154 INTEGER,
PARAMETER:: SOLWEIG_ldown = 0
161 ALLOCATE (knight(1, 1))
162 ALLOCATE (tgmap0(1, 1))
163 ALLOCATE (svf_blgd_veg(1, 1))
169 IF (posture == 1)
THEN 184 height = heightgravity
222 tmp = 1 - (svf + svfveg - 1)
223 WHERE (tmp <= 0) tmp = 0.000000001
224 svfalfa = asin(exp(log(tmp)/2))
227 first = int(anint(height))
231 second = int(anint(height*20))
234 svf_blgd_veg = (svf - (1 - svfveg)*(1 - psi))
237 CALL daylen(doy, lat, dayl, dec, sndn, snup)
239 altitude = 90 - zenith_deg
242 ea = 6.107*10**((7.5*ta)/(237.3 + ta))*(rh/100)
243 msteg = 46.5*(ea/(ta + 273.15))
244 emis_sky = (1 - (1 + msteg)*exp(-((1.2 + 3.0*msteg)**0.5))) - 0.04
247 IF (altitude > 0)
THEN 251 CALL clearnessindex_2013b(zen, doy, ta, rh/100, radg, lat, p, i0, ci, kt, i0et, ciuncorr)
256 IF (onlyglobal == 1)
THEN 261 IF (usevegdem == 1)
THEN 273 DO j = 1,
SIZE(azimutha)
274 CALL sunonsurface_veg(azimutha(j), scale, buildings, first, second, psi, sos)
277 gvf = gvf/
SIZE(azimutha) + (buildings*(-1) + 1)
284 kdown2d = radi*shadow*sin(altitude*
deg2rad) &
285 + radd*svf_blgd_veg &
286 + radg*alb_bldg*(1 - svf_blgd_veg)*(1 - f_sh)
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))
295 radi, radg, radd, azimuth, altitude, psi, t, alb_ground, &
296 keast, knorth, ksouth, kwest)
336 timeadd = (dectime + 1) - doy - (it/24.)
337 weight1 = exp(-33.27*timeadd)
338 tgmap1 = tgmap0*(1 - weight1) + tgmap1*weight1
339 lup2d = sbc*emis_ground*((tgmap1 + 273.15)**4)
345 IF (dectime < (doy + 0.5) .AND. dectime > doy .AND. altitude < 1.0)
THEN 384 lup2d = sbc*emis_ground*((knight + ta + tg + 273.15)**4)
389 IF (solweig_ldown == 1)
THEN 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)
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))
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
414 altitude, ta, tw, sbc, emis_wall, emis_sky, t, ci, azimuth, zen, ldown, svfalfa, &
415 least, lnorth, lsouth, lwest)
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
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]
441 DEALLOCATE (svf_blgd_veg)
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
466 REAL(KIND(1d0)),
INTENT(in)::lamdaP
467 REAL(KIND(1d0)),
INTENT(in)::lamdaF
469 REAL(KIND(1d0))::lamdaW
473 hw = (lamdaw*lamdap)/(2*lamdap*(1 - lamdap))
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
484 REAL(KIND(1d0)),
INTENT(in)::hw
485 REAL(KIND(1d0))::svfGround
488 svfground = a*exp(b*hw) + c*exp(d*hw)
494 REAL(KIND(1d0)),
PARAMETER::e = 0.5572
495 REAL(KIND(1d0)),
PARAMETER::f = 0.0589
496 REAL(KIND(1d0)),
PARAMETER::g = 0.4143
498 REAL(KIND(1d0)),
INTENT(in)::hw
499 REAL(KIND(1d0))::svfRoof
502 svfroof = e*exp(-f*hw) + g
507 SUBROUTINE clearnessindex_2013b(zen, DOY, Ta, RH, radG, lat, P_kPa, I0, CI, Kt, I0et, CIuncorr)
513 INTEGER,
intent(in):: DOY
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
546 IF (p_kpa == -999)
THEN 556 m = 35.*cos(zen)*((1224.*(cos(zen)**2) + 1.)**(-1./2.))
557 trpg = 1.021 - 0.084*(m*(0.000949*p + 0.051))**0.5
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]
580 IF (doy > 335 .OR. doy <= 60)
THEN 582 ELSE IF (doy > 60 .AND. doy <= 152)
THEN 584 ELSE IF (doy > 152 .AND. doy <= 244)
THEN 586 ELSE IF (doy > 244 .AND. doy <= 335)
THEN 592 td = (b2*(((a2*ta)/(b2 + ta)) + log(rh)))/(a2 - (((a2*ta)/(b2 + ta)) + log(rh)))
594 u = exp(0.1133 - log(ig + 1) + 0.0393*td)
595 tw = 1 - 0.077*((u*m)**0.3)
598 i0 = itoa*cos(zen)*trpg*tw*d*tar
608 corr = 0.1473*log(90 - (zen/pi*180)) + 0.3454
611 ci = ciuncorr + (1 - corr)
612 i0et = itoa*cos(zen)*d
630 REAL(KIND(1d0)) ::b, D
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))
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
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
656 ALLOCATE (alfa(1, 1))
663 ALLOCATE (ukil(1, 1))
665 ALLOCATE (ssurf(1, 1))
666 ALLOCATE (hkil(1, 1))
671 xa = 1.-2./(tan(alfa)*tan(beta))
672 ha = 2./(tan(alfa)*tan(beta))
691 za = (ba**2 - qa**2/4.)**0.5
695 ai = (sin(phi) - phi*cos(phi))/(1 - cos(phi))
702 f_sh = (2*pi*ba - ssurf)/(2*pi*ba)
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
729 REAL(KIND(1d0)),
intent(out)::radI
730 REAL(KIND(1d0))::alfa
735 IF (ta <= -99 .OR. rh <= -99)
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 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))
754 radd = max(0.d0, radd)
755 radd = min(radg, radd)
758 radi = (radg - radd)/(sin(alfa))
768 IF (altitude < 1 .AND. radi > radg)
THEN 1047 radI, radG, radD, azimuth, altitude, psi, t, albedo, & ! input
1048 Keast, Knorth, Ksouth, Kwest)
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
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
1065 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW
1066 REAL(KIND(1d0)),
DIMENSION(1, 1) :: viktveg, viktwall
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
1079 REAL(KIND(1d0)),
ALLOCATABLE,
DIMENSION(:, :) :: svfviktbuveg
1081 ALLOCATE (svfviktbuveg(1, 1))
1085 azis = azimuth - 90 + t
1086 aziw = azimuth - 180 + t
1087 azin = azimuth - 270 + t
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)
1096 keast = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh)
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)
1105 ksouth = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh)
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)
1114 kwest = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh)
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)
1123 knorth = radd*(1 - svfviktbuveg) + radg*albedo*svfviktbuveg*(1 - f_sh)
1126 DEALLOCATE (svfviktbuveg)
1129 SUBROUTINE kvikt_veg(isvf, isvfveg, vikttot, & !input
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
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
1145 svfvegbu = (isvfveg + isvf - 1)
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
1155 altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, zen, ldown, svfalfa, &
1156 Least, Lnorth, Lsouth, Lwest)
1160 REAL(KIND(1D0)),
intent(in)::altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, ldown, zen
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
1166 REAL(KIND(1D0))::vikttot, aziE, aziN, aziS, aziW, c
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
1176 REAL(KIND(1d0)),
DIMENSION(1, 1) :: viktveg, viktsky, viktrefl, viktwall
1177 REAL(KIND(1d0)),
DIMENSION(1, 1) :: F_sh
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
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))
1221 oneminussvfe = 1.-svfe;
WHERE (oneminussvfe <= 0) oneminussvfe = 0.000000001
1222 oneminussvfs = 1.-svfs;
WHERE (oneminussvfs <= 0) oneminussvfs = 0.000000001
1223 oneminussvfw = 1.-svfw;
WHERE (oneminussvfw <= 0) oneminussvfw = 0.000000001
1224 oneminussvfn = 1.-svfn;
WHERE (oneminussvfn <= 0) oneminussvfn = 0.000000001
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))
1232 azin = azimuth - 90 + t
1233 azie = azimuth - 180 + t
1234 azis = azimuth - 270 + t
1240 IF (solweig_ldown == 1)
THEN 1242 lsky_allsky = emis_sky*sbc*((ta + 273.15)**4)*(1 - c) + c*sbc*((ta + 273.15)**4)
1248 CALL lvikt_veg(svfe, svfeveg, svfeaveg, vikttot, &
1249 viktveg, viktsky, viktrefl, viktwall)
1251 IF (altitude > 0)
THEN 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
1261 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1265 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1267 lsky = ((svfe + svfeveg - 1)*lsky_allsky)*viktsky*0.5
1268 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1270 lrefl = (ldown2d+lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1271 least = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1274 CALL lvikt_veg(svfs, svfsveg, svfsaveg, vikttot, &
1275 viktveg, viktsky, viktrefl, viktwall)
1277 IF (altitude > 0)
THEN 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
1287 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1291 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1293 lsky = ((svfs + svfsveg - 1)*lsky_allsky)*viktsky*0.5
1294 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1296 lrefl = (ldown2d+lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1297 lsouth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1300 CALL lvikt_veg(svfw, svfwveg, svfwaveg, vikttot, &
1301 viktveg, viktsky, viktrefl, viktwall)
1303 IF (altitude > 0)
THEN 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
1313 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1317 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1319 lsky = ((svfw + svfwveg - 1)*lsky_allsky)*viktsky*0.5
1320 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1322 lrefl = (ldown2d+lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1323 lwest = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1326 CALL lvikt_veg(svfn, svfnveg, svfnaveg, vikttot, &
1327 viktveg, viktsky, viktrefl, viktwall)
1329 IF (altitude > 0)
THEN 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
1339 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1343 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1345 lsky = ((svfn + svfnveg - 1)*lsky_allsky)*viktsky*0.5
1346 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1348 lrefl = (ldown2d+lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1349 lnorth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1351 DEALLOCATE (svfalfae)
1352 DEALLOCATE (svfalfas)
1353 DEALLOCATE (svfalfaw)
1354 DEALLOCATE (svfalfan)
1357 DEALLOCATE (betasun)
1358 DEALLOCATE (lground)
1361 DEALLOCATE (lsky_allsky)
1363 DEALLOCATE (lwallsh)
1364 DEALLOCATE (lwallsun)
1365 DEALLOCATE (viktonlywall)
1366 DEALLOCATE (viktaveg)
1367 DEALLOCATE (svfvegbu)
1371 SUBROUTINE lvikt_veg(isvf, isvfveg, isvfaveg, vikttot, & ! input
1372 viktveg, viktsky, viktrefl, viktwall)
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
1384 REAL(KIND(1d0)),
DIMENSION(1, 1) :: viktonlywall
1385 REAL(KIND(1d0)),
DIMENSION(1, 1) :: viktaveg
1386 REAL(KIND(1d0)),
DIMENSION(1, 1) :: svfvegbu
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
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
1405 viktwall = viktonlywall - viktaveg
1407 svfvegbu = (isvfveg + isvf - 1)
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
1424 SUBROUTINE issign(IX, MAXPOS, ISIGNM)
1425 REAL(KIND(1d0)) IX, MAXPOS, ISIGNM
1427 IF (ix < 0 .OR. ix > maxpos) isignm = -1
1428 IF (ix == 0) isignm = 0
1434 CHARACTER(len=20) FUNCTION str(k)
1435 INTEGER,
INTENT(in) :: k
1530 REAL(KIND(1d0)),
intent(in) :: azimuth, altitude, scale
1531 REAL(KIND(1d0)),
DIMENSION(1, 1),
intent(out) :: shadow
1533 REAL(KIND(1d0)),
DIMENSION(1, 1):: a
1535 REAL(KIND(1d0)),
PARAMETER :: pi = 3.141592653589793
1536 REAL(KIND(1d0)),
PARAMETER :: maxpos = 10000000000.0
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
1544 REAL(KIND(1d0)),
ALLOCATABLE,
DIMENSION(:, :) :: temp, f
1559 azi = min(azimuth, 0 - 0.0001)*degrees
1560 alt = min(altitude, 90 - 0.0001)*degrees
1563 ALLOCATE (temp(1, 1))
1577 amaxvalue = maxval(a)
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)
1589 dssin = abs(1./sinazimuth)
1590 dscos = abs(1./cosazimuth)
1591 tanaltitudebyscale = tan(alt)/scale
1594 DO WHILE (amaxvalue >= dz .AND. abs(dx) <= 1 .AND. abs(dy) <= 1)
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))
1601 dy = signsinazimuth*abs(nint(index*tanazimuth))
1602 dx = -1.*signcosazimuth*index
1606 dz = ds*index*tanaltitudebyscale
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))
1621 temp(xp1:xp2, yp1:yp2) = a(xc1:xc2, yc1:yc2) - dz
1860 SUBROUTINE sunonsurface_veg(iazimuthA, scale, buildings, first, second, psi, sos)
1868 INTEGER::first, second
1870 REAL(KIND(1d0)),
intent(in)::iazimuthA
1871 REAL(KIND(1d0)),
intent(in)::scale
1872 REAL(KIND(1d0)),
intent(in)::psi
1874 REAL(KIND(1d0)),
DIMENSION(1, 1),
intent(in) :: buildings
1875 REAL(KIND(1d0)),
DIMENSION(1, 1),
intent(out)::sos
1877 REAL(KIND(1d0)),
DIMENSION(1, 1):: sunwall
1879 REAL(KIND(1d0)),
DIMENSION(1, 1):: sh, vegsh
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
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
1890 REAL(KIND(1d0)),
PARAMETER :: pi = 3.141592653589793
1891 REAL(KIND(1d0)),
PARAMETER :: maxpos = 10000000000.0
1893 ALLOCATE (weightsumwall(1, 1))
1894 ALLOCATE (weightsumsh(1, 1))
1895 ALLOCATE (gvf1(1, 1))
1896 ALLOCATE (gvf2(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))
1905 iazimuth = iazimutha*(pi/180)
1907 IF (iazimuth == 0)
THEN 1908 iazimuth = iazimuth + 0.000001
1919 sh1 = sh - (1 - vegsh)*(1 - psi)
1930 weightsumwall = 0.0d0
1932 first = int(
REAL(first, kind(1d0))*scale)
1933 second = int(
REAL(second, kind(1d0))*scale)
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)
1947 dssin = abs(1./sinazimuth)
1948 dscos = abs(1./cosazimuth)
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))
1958 dy = signsinazimuth*abs(nint(index*tanazimuth))
1959 dx = -1.*signcosazimuth*index
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))
1975 tempbu(xp1:xp2, yp1:yp2) = buildings(xc1:xc2, yc1:yc2)
1977 tempsh(xp1:xp2, yp1:yp2) = sh1(xc1:xc2, yc1:yc2)
1980 weightsumsh = weightsumsh + tempsh*f
1982 tempwallsun(xp1:xp2, yp1:yp2) = sunwall(xc1:xc2, yc1:yc2)
1983 tempb = tempwallsun*f
1984 WHERE ((tempb + tempbub) > 0)
1988 weightsumwall = weightsumwall + tempbub
1990 IF (index*scale == first)
THEN 1991 gvf1 = (weightsumwall + weightsumsh)/first
1999 gvf2 = (weightsumsh + weightsumwall)/second
2005 sos = (gvf1*0.5 + gvf2*0.4)/0.9
2007 DEALLOCATE (weightsumwall)
2008 DEALLOCATE (weightsumsh)
2014 DEALLOCATE (tempbub)
2015 DEALLOCATE (tempwallsun)
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)
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)