31 REAL(kind(1d0)),
PARAMETER ::
pi = atan(1.)*4
37 SUBROUTINE beers_cal_main(iy, id, dectime, lamdaP, lamdaF, avkdn, ldown, Temp_C, avrh, &
38 Press_hPa, Tsurf, lat, lng, alt, timezone, zenith_deg, azimuth, &
39 alb_ground, alb_bldg, emis_ground, emis_wall, &
45 INTEGER,
INTENT(in) :: iy
46 INTEGER,
INTENT(in) :: id
47 REAL(KIND(1D0)),
INTENT(in) :: lamdaP
48 REAL(KIND(1D0)),
INTENT(in) :: lamdaF
53 REAL(KIND(1D0)),
INTENT(in) :: Press_hPa
54 REAL(KIND(1D0)),
INTENT(in) :: Temp_C
55 REAL(KIND(1D0)),
INTENT(in) :: avrh
56 REAL(KIND(1D0)),
INTENT(in) :: avkdn
57 REAL(KIND(1D0)),
INTENT(in) :: ldown
58 REAL(KIND(1D0)),
INTENT(in) :: Tsurf
61 REAL(KIND(1D0)),
INTENT(in) :: zenith_deg
62 REAL(KIND(1D0)),
INTENT(in) :: azimuth
63 REAL(KIND(1D0)),
INTENT(in) :: dectime
64 REAL(KIND(1D0)),
INTENT(in) :: timezone, lat, lng, alt
65 REAL(KIND(1D0)),
INTENT(in) :: alb_bldg
66 REAL(KIND(1D0)),
INTENT(in) :: alb_ground
67 REAL(KIND(1D0)),
INTENT(in) :: emis_wall
68 REAL(KIND(1D0)),
INTENT(in) :: emis_ground
70 REAL(KIND(1D0)),
PARAMETER :: absL = 0.97
71 REAL(KIND(1D0)),
PARAMETER :: absK = 0.7
72 REAL(KIND(1D0)),
PARAMETER :: Fside = 0.22
73 REAL(KIND(1D0)),
PARAMETER :: Fup = 0.06
76 REAL(KIND(1D0)),
DIMENSION(ncolumnsDataOutBEERS - 5),
INTENT(OUT) :: dataOutLineBEERS
78 REAL(KIND(1D0)) :: t, psi
79 REAL(KIND(1D0)) :: altitude, zen
80 REAL(KIND(1D0)) :: CI, c, I0, Kt, Tw, Tg
81 REAL(KIND(1D0)) :: Ta, RH, P, radG, radD, radI
82 REAL(KIND(1D0)) :: I0et, CIuncorr
83 REAL(KIND(1D0)) :: SNDN, SNUP, DEC, DAYL
84 REAL(KIND(1D0)) :: msteg, emis_sky, ea
85 REAL(KIND(1D0)) :: shadowground, shadowwalls, shadowroof
88 REAL(KIND(1D0)) :: CIlatenight
89 REAL(KIND(1D0)) :: dectime_sunrise, zen_sunrise, I0_sunrise
94 REAL(KIND(1D0)) :: svfalfa
96 REAL(KIND(1D0)) :: Tmrt, Sstr, F_sh
98 REAL(KIND(1D0)) :: tmp, altmax
99 REAL(KIND(1D0)) :: svf_bldg_veg
100 REAL(KIND(1D0)) :: svf_ground, svf_roof
101 REAL(KIND(1D0)) :: svf_veg
102 REAL(KIND(1D0)) :: svf_aveg
103 REAL(KIND(1D0)) :: Kdown, Keast, Knorth, Ksouth, Kup2d, Kwest
104 REAL(KIND(1D0)) :: Ldown2d, Least, Lnorth, Lsouth, Lup2d, Lwest
111 REAL(KIND(1D0)),
PARAMETER :: SBC = 5.67051e-8
113 INTEGER,
PARAMETER :: onlyglobal = 1
120 INTEGER,
PARAMETER :: SOLWEIG_ldown = 1
121 INTEGER,
PARAMETER :: BEERS_tsurf = 1
165 tmp = 1 - (svf_ground + svf_veg - 1)
166 IF (tmp <= 1.e-6) tmp = 1.e-6
167 svfalfa = asin(exp(log(tmp)/2))
170 svf_bldg_veg = (svf_ground - (1 - svf_veg)*(1 - psi))
173 CALL daylen(doy, lat, dayl, dec, sndn, snup)
175 altitude = 90 - zenith_deg
178 ea = 6.107*10**((7.5*ta)/(237.3 + ta))*(rh/100)
179 msteg = 46.5*(ea/(ta + 273.15))
180 emis_sky = (1 - (1 + msteg)*exp(-((1.2 + 3.0*msteg)**0.5)))
183 IF (altitude > 0.1)
THEN
188 i0, ci, kt, i0et, ciuncorr)
192 IF (onlyglobal == 1)
THEN
202 CALL kroof(radi, radd, radg, f_sh, altitude, svf_roof, svf_veg, shadowroof, psi, alb_bldg, kdown)
208 kup2d = alb_ground*( &
209 radi*shadowground*sin(altitude*
deg2rad) &
210 + radd*svf_bldg_veg &
211 + alb_bldg*(1 - svf_bldg_veg)*(radg*(1 - f_sh) + radd*f_sh))
215 svf_ground, svf_veg, shadowground, f_sh, &
216 radi, radg, radd, azimuth, altitude, psi, t, alb_ground, alb_bldg, &
217 keast, knorth, ksouth, kwest)
219 IF (beers_tsurf == 1)
THEN
220 CALL tsurfbeers(iy, ta, rh, radi, i0, dectime, snup, altitude, zen, timezone, lat, lng, alt, &
229 lup2d = sbc*emis_ground*((shadowground*tg + ta + 273.15)**4)
234 cilatenight, dectime_sunrise, zen_sunrise, i0_sunrise)
242 IF (beers_tsurf == 1)
THEN
254 IF (dectime < (doy + 0.5) .AND. dectime > doy .AND. altitude < 1.0)
THEN
288 lup2d = sbc*emis_ground*((ta + tg + 273.15)**4)
293 IF (solweig_ldown == 1)
THEN
294 ldown2d = (svf_roof + svf_veg - 1)*emis_sky*sbc*((ta + 273.15)**4) &
295 + (2 - svf_veg - svf_aveg)*emis_wall*sbc*((ta + 273.15)**4) &
296 + (svf_aveg - svf_roof)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
297 + (2 - svf_roof - svf_veg)*(1 - emis_wall)*emis_sky*sbc*((ta + 273.15)**4)
301 ldown2d = ldown2d*(1 - c) + c*( &
302 (svf_roof + svf_veg - 1)*sbc*((ta + 273.15)**4) &
303 + (2 - svf_veg - svf_aveg)*emis_wall*sbc*((ta + 273.15)**4) &
304 + (svf_aveg - svf_roof)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
305 + (2 - svf_roof - svf_veg)*(1 - emis_wall)*sbc*((ta + 273.15)**4))
309 ldown2d = (svf_roof + svf_veg - 1)*ldown &
310 + (2 - svf_veg - svf_aveg)*emis_wall*sbc*((ta + 273.15)**4) &
311 + (svf_aveg - svf_roof)*emis_wall*sbc*((ta + 273.15 + tw)**4) &
312 + (2 - svf_roof - svf_veg)*(1 - emis_wall)*ldown
316 CALL lwalls(svf_ground, svf_veg, svf_aveg, &
318 altitude, ta, tw, sbc, emis_wall, &
319 emis_sky, t, ci, azimuth, ldown, svfalfa, f_sh, &
320 least, lnorth, lsouth, lwest)
323 sstr = absk*(kdown*fup + kup2d*fup + knorth*fside + keast*fside + ksouth*fside + kwest*fside) &
324 + absl*(ldown2d*fup + lup2d*fup + lnorth*fside + least*fside + lsouth*fside + lwest*fside)
325 tmrt = sqrt(sqrt((sstr/(absl*sbc)))) - 273.15
327 dataoutlinebeers = [azimuth, altitude, radg, radi, radd, &
328 kdown, kup2d, ksouth, kwest, knorth, keast, &
329 ldown2d, lup2d, lsouth, lwest, lnorth, least, &
330 tmrt, i0, ci, shadowground, shadowwalls, svf_ground, svf_roof, svf_bldg_veg, &
340 CIlatenight, dectime_sunrise, zen_sunrise, I0_sunrise)
346 INTEGER,
INTENT(in) :: iy
347 INTEGER,
INTENT(in) :: DOY
349 REAL(KIND(1D0)),
INTENT(in) :: Ta_degC
350 REAL(KIND(1D0)),
INTENT(in) :: RH_frac
351 REAL(KIND(1D0)),
INTENT(in) :: P_kPa
352 REAL(KIND(1D0)),
INTENT(in) :: radG
353 REAL(KIND(1D0)),
INTENT(in) :: lat
354 REAL(KIND(1D0)),
INTENT(out) :: CIlatenight
355 REAL(KIND(1D0)),
INTENT(out) :: dectime_sunrise
356 REAL(KIND(1D0)),
INTENT(out) :: zen_sunrise
357 REAL(KIND(1D0)),
INTENT(out) :: I0_sunrise
358 REAL(KIND(1D0)) :: I0et
359 REAL(KIND(1D0)) :: CIuncorr
360 REAL(KIND(1D0)) :: Kt
361 REAL(KIND(1D0)) :: DAYL, DEC, SNDN, SNUP, azimuth
363 CALL daylen(doy, lat, dayl, dec, sndn, snup)
364 dectime_sunrise = doy + (snup + .5)/24.
366 REAL(iy, KIND(1D0)), &
368 0.D0, lat, 0.D0, 100.D0, &
369 azimuth, zen_sunrise)
371 zen_sunrise = zen_sunrise/180.*
pi
374 i0_sunrise, cilatenight, kt, i0et, ciuncorr)
379 radI, radD, radG, F_sh, altitude, svfr, svfveg, shadow, psi, alb_bldg, & ! input
384 REAL(KIND(1D0)),
INTENT(in) :: radI
385 REAL(KIND(1D0)),
INTENT(in) :: radG
386 REAL(KIND(1D0)),
INTENT(in) :: radD
388 REAL(KIND(1D0)),
INTENT(in) :: altitude
389 REAL(KIND(1D0)),
INTENT(in) :: svfr
390 REAL(KIND(1D0)),
INTENT(in) :: psi
391 REAL(KIND(1D0)),
INTENT(in) :: alb_bldg
392 REAL(KIND(1D0)),
INTENT(in) :: shadow, svfveg
393 REAL(KIND(1D0)),
INTENT(in) :: F_sh
394 REAL(KIND(1D0)),
INTENT(out) :: Kdown
397 REAL(KIND(1D0)) :: svfrbuveg
414 svfrbuveg = svfr - (1.0 - svfveg)*(1.0 - psi)
416 kdown = radi*shadow*sin(altitude*
deg2rad) &
418 + alb_bldg*(1 - svfrbuveg)*(radg*(1 - f_sh) + radd*f_sh)
425 REAL(KIND(1D0)),
INTENT(in) :: svfr
426 REAL(KIND(1D0)),
INTENT(in) :: svfveg
427 REAL(KIND(1D0)),
INTENT(out) :: svfalfa
428 REAL(KIND(1D0)),
INTENT(out) :: tmp
432 tmp = 1.-(svfr + svfveg - 1.)
433 IF (tmp <= 0) tmp = 0.000000001
434 svfalfa = asin(exp(log(tmp)/2.))
440 REAL(kind(1d0)),
PARAMETER :: a = 0.5598
441 REAL(kind(1d0)),
PARAMETER :: b = -0.2485
442 REAL(kind(1d0)),
PARAMETER :: c = 0.4112
443 REAL(kind(1d0)),
PARAMETER :: d = -0.02528
445 REAL(kind(1d0)),
INTENT(in) :: lamdap
446 REAL(kind(1d0)),
INTENT(in) :: lamdaf
447 REAL(kind(1d0)) :: hw
448 REAL(kind(1d0)) :: lamdaw
453 hw = (lamdaw*lamdap)/(2*lamdap*(1 - lamdap))
459 REAL(kind(1d0)),
PARAMETER :: a = 0.5598
460 REAL(kind(1d0)),
PARAMETER :: b = -0.2485
461 REAL(kind(1d0)),
PARAMETER :: c = 0.4112
462 REAL(kind(1d0)),
PARAMETER :: d = -0.02528
464 REAL(kind(1d0)),
INTENT(in) :: hw
465 REAL(kind(1d0)) :: svfground
468 svfground = a*exp(b*hw) + c*exp(d*hw)
474 REAL(kind(1d0)),
PARAMETER :: e = 0.5572
475 REAL(kind(1d0)),
PARAMETER :: f = 0.0589
476 REAL(kind(1d0)),
PARAMETER :: g = 0.4143
478 REAL(kind(1d0)),
INTENT(in) :: hw
479 REAL(kind(1d0)) :: svfroof
482 svfroof = e*exp(-f*hw) + g
486 SUBROUTINE tsurfbeers(iy, Ta, RH, radI, I0, dectime, SNUP, altitude, zen, timezone, lat, lng, alt, &
492 INTEGER,
INTENT(in) :: iy
493 REAL(KIND(1D0)),
INTENT(in) :: Ta, RH, radI, I0
494 REAL(KIND(1D0)),
INTENT(in) :: dectime, SNUP, altitude, zen
495 REAL(KIND(1D0)),
INTENT(in) :: timezone, lat, lng, alt
496 REAL(KIND(1D0)),
INTENT(out) :: Tg, Tgwall, altmax
497 REAL(KIND(1D0)),
PARAMETER :: TgK = 0.37
498 REAL(KIND(1D0)),
PARAMETER :: Tstart = 3.41
499 REAL(KIND(1D0)),
PARAMETER :: TgK_wall = 0.37
500 REAL(KIND(1D0)),
PARAMETER :: Tstart_wall = -3.41
501 REAL(KIND(1D0)),
PARAMETER :: TmaxLST = 15.
502 REAL(KIND(1D0)),
PARAMETER :: TmaxLST_wall = 15.
504 REAL(KIND(1D0)) :: Ktc, notU, Tgamp, Tgampwall, radI0, corr, CI_Tg
505 REAL(KIND(1D0)) :: fifteen, sunmaximum, zen_sunmax, dectimemax, azimuth
511 DO WHILE (sunmaximum <= 90.-zen_sunmax)
512 sunmaximum = 90.-zen_sunmax
513 fifteen = fifteen + 15./1440.
514 dectimemax = floor(dectime) + 10/24.+fifteen
516 REAL(iy, KIND(1D0)), &
518 timezone, lat, lng, alt, &
524 tgamp = (tgk*altmax - tstart) + tstart
525 tgampwall = (tgk_wall*altmax - (tstart_wall)) + (tstart_wall)
526 tg = tgamp*sin((((dectime - floor(dectime)) - snup/24)/(tmaxlst/24 - snup/24))*
pi/2) + tstart
527 tgwall = tgampwall*sin((((dectime - floor(dectime)) - snup/24)/(tmaxlst_wall/24 - snup/24))*
pi/2) + (tstart_wall)
531 corr = 0.1473*log(90.-(zen*
rad2deg)) + 0.3454
532 ci_tg = (radi/radi0) + (1.-corr)
535 IF (ci_tg > 1) ci_tg = 1
537 tgwall = tgwall*ci_tg
539 IF (tgwall < 0) tgwall = 0
552 REAL(KIND(1D0)),
INTENT(in) :: HW, azimuth, zen
553 REAL(KIND(1D0)),
INTENT(out) :: shadowground, shadowwalls
555 REAL(KIND(1D0)) :: ROOF_HEIGHT, ROAD_WIDTH, THETA_Z, THETA_S, THETA_S180, ndir, wx
556 REAL(KIND(1D0)) :: LshadowRoad(1:8), Wallsun(1:8)
563 theta_s180 = abs(
pi - theta_s)
569 wx = 1./sin(abs((float(idir)*
pi/ndir) - theta_s180))
570 lshadowroad(idir) = (roof_height/road_width)*tan(theta_z)*sin(abs((float(idir)*
pi/ndir) - theta_s180))
571 wallsun(idir) = (wx/tan(theta_z))/(roof_height/road_width)
575 WHERE (lshadowroad >= 1.) lshadowroad = 1.
576 lshadowroad = -1*lshadowroad + 1.
577 shadowground = sum(lshadowroad)/ndir
580 WHERE (wallsun >= 1.) wallsun = 1.
581 shadowwalls = (sum(wallsun)/ndir)/2.
586 I0, CI, Kt, I0et, CIuncorr)
590 INTEGER,
INTENT(in) :: DOY
591 REAL(KIND(1D0)),
INTENT(in) :: zen
592 REAL(KIND(1D0)),
INTENT(in) :: Ta_degC
593 REAL(KIND(1D0)),
INTENT(in) :: RH_frac
594 REAL(KIND(1D0)),
INTENT(in) :: P_kPa
595 REAL(KIND(1D0)),
INTENT(in) :: radG
596 REAL(KIND(1D0)),
INTENT(in) :: lat
597 REAL(KIND(1D0)),
INTENT(out) :: I0
598 REAL(KIND(1D0)),
INTENT(out) :: CI
599 REAL(KIND(1D0)),
INTENT(out) :: Kt
600 REAL(KIND(1D0)),
INTENT(out) :: I0et
601 REAL(KIND(1D0)),
INTENT(out) :: CIuncorr
602 REAL(KIND(1D0)) :: iG, Itoa, p
603 REAL(KIND(1D0)),
DIMENSION(4) :: G
619 IF (p_kpa == -999)
THEN
627 m = 35.*cos(zen)*((1224.*(cos(zen)**2) + 1.)**(-1./2.))
628 trpg = 1.021 - 0.084*(m*(0.000949*p + 0.051))**0.5
633 g = [3.37, 2.85, 2.80, 2.64]
634 ELSE IF (lat >= 10 .AND. lat < 20)
THEN
635 g = [2.99, 3.02, 2.70, 2.93]
636 ELSE IF (lat >= 20 .AND. lat < 30)
THEN
637 g = [3.60, 3.00, 2.98, 2.93]
638 ELSE IF (lat >= 30 .AND. lat < 40)
THEN
639 g = [3.04, 3.11, 2.92, 2.94]
640 ELSE IF (lat >= 40 .AND. lat < 50)
THEN
641 g = [2.70, 2.95, 2.77, 2.71]
642 ELSE IF (lat >= 50 .AND. lat < 60)
THEN
643 g = [2.52, 3.07, 2.67, 2.93]
644 ELSE IF (lat >= 60 .AND. lat < 70)
THEN
645 g = [1.76, 2.69, 2.61, 2.61]
646 ELSE IF (lat >= 70 .AND. lat < 80)
THEN
647 g = [1.60, 1.67, 2.24, 2.63]
648 ELSE IF (lat >= 80 .AND. lat < 90)
THEN
649 g = [1.11, 1.44, 1.94, 2.02]
651 IF (doy > 335 .OR. doy <= 60)
THEN
653 ELSE IF (doy > 60 .AND. doy <= 152)
THEN
655 ELSE IF (doy > 152 .AND. doy <= 244)
THEN
657 ELSE IF (doy > 244 .AND. doy <= 335)
THEN
663 td = (b2*(((
a2*ta_degc)/(b2 + ta_degc)) + log(rh_frac)))/(
a2 - (((
a2*ta_degc)/(b2 + ta_degc)) + log(rh_frac)))
665 u = exp(0.1133 - log(ig + 1) + 0.0393*td)
666 tw = 1 - 0.077*((u*m)**0.3)
669 i0 = itoa*cos(zen)*trpg*tw*d*tar
679 corr = 0.1473*log(90 - (zen/
pi*180)) + 0.3454
682 ci = ciuncorr + (1 - corr)
683 i0et = itoa*cos(zen)*d
699 REAL(KIND(1D0)) :: b, D
701 b = 2*3.141592654*jday/365
702 d = sqrt(1.00011 + 0.034221*cos(b) + 0.001280*sin(b) + 0.000719*cos(2*b) + 0.000077*sin(2*b))
712 REAL(KIND(1D0)),
INTENT(in) :: zen
713 REAL(KIND(1D0)),
INTENT(in) :: svfalfa
714 REAL(KIND(1D0)),
INTENT(out) :: F_sh
715 REAL(KIND(1D0)) :: beta
716 REAL(KIND(1D0)) :: alfa, xa, ha, hkil, ba
717 REAL(KIND(1D0)) :: Ai, phi, qa, Za
718 REAL(KIND(1D0)) :: ukil, Ssurf
735 xa = 1.-2./(tan(alfa)*tan(beta))
736 ha = 2./(tan(alfa)*tan(beta))
741 IF (xa < 0) qa = tan(beta)/2
748 za = (ba**2 - qa**2/4.)**0.5
750 ai = (sin(phi) - phi*cos(phi))/(1 - cos(phi))
756 f_sh = (2*
pi*ba - ssurf)/(2*
pi*ba)
758 IF (f_sh < 0) f_sh = 0.0
759 IF (f_sh > 0.5) f_sh = 0.5
780 REAL(KIND(1D0)),
INTENT(in) :: radG
781 REAL(KIND(1D0)),
INTENT(in) :: altitude
782 REAL(KIND(1D0)),
INTENT(in) :: Kt
783 REAL(KIND(1D0)),
INTENT(in) :: Ta
784 REAL(KIND(1D0)),
INTENT(in) :: RH
785 REAL(KIND(1D0)),
INTENT(out) :: radD
786 REAL(KIND(1D0)),
INTENT(out) :: radI
787 REAL(KIND(1D0)) :: alfa
791 IF (ta <= -99 .OR. rh <= -99)
THEN
793 radd = radg*(1.020 - 0.248*kt)
794 ELSE IF (kt > 0.3 .AND. kt < 0.78)
THEN
795 radd = radg*(1.45 - 1.67*kt)
796 ELSE IF (kt >= 0.78)
THEN
801 radd = radg*(1 - 0.232*kt + 0.0239*sin(alfa) - 0.000682*ta + 0.0195*(rh/100))
802 ELSE IF (kt > 0.3 .AND. kt < 0.78)
THEN
803 radd = radg*(1.329 - 1.716*kt + 0.267*sin(alfa) - 0.00357*ta + 0.106*(rh/100))
804 ELSE IF (kt >= 0.78)
THEN
805 radd = radg*(0.426*kt - 0.256*sin(alfa) + 0.00349*ta + 0.0734*(rh/100))
809 radd = max(0.d0, radd)
810 radd = min(radg, radd)
813 radi = (radg - radd)/(sin(alfa))
816 IF (altitude < 1 .AND. radi > radg)
THEN
823 svf, svfveg, shadow, F_sh, &
824 radI, radG, radD, azimuth, altitude, psi, t, alb_ground, alb_bldg, & ! input
825 Keast, Knorth, Ksouth, Kwest)
830 REAL(KIND(1D0)),
INTENT(in) :: radI
831 REAL(KIND(1D0)),
INTENT(in) :: radG
832 REAL(KIND(1D0)),
INTENT(in) :: radD
833 REAL(KIND(1D0)),
INTENT(in) :: azimuth
834 REAL(KIND(1D0)),
INTENT(in) :: altitude
835 REAL(KIND(1D0)),
INTENT(in) :: psi
836 REAL(KIND(1D0)),
INTENT(in) :: t
837 REAL(KIND(1D0)),
INTENT(in) :: alb_bldg
838 REAL(KIND(1D0)),
INTENT(in) :: alb_ground
841 REAL(KIND(1D0)),
INTENT(in) :: shadow, F_sh, svf, svfveg
842 REAL(KIND(1D0)),
INTENT(out) :: Keast, Knorth, Ksouth, Kwest
844 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW
845 REAL(KIND(1D0)) :: viktveg, viktwall
846 REAL(KIND(1D0)) :: KeastI, KsouthI, KwestI, KnorthI, Kuptowall
847 REAL(KIND(1D0)) :: KeastDG, KsouthDG, KwestDG, KnorthDG
848 REAL(KIND(1D0)) :: svfE, svfS, svfW, svfN, svfEveg, svfSveg, svfWveg, svfNveg, svfbuveg
849 REAL(KIND(1D0)) :: gvfalb, gvfalbnosh
852 REAL(KIND(1D0)) :: svfviktbuveg
857 azis = azimuth - 90 + t
858 aziw = azimuth - 180 + t
859 azin = azimuth - 270 + t
871 IF (azimuth > (360 - t) .OR. azimuth <= (180 - t))
THEN
876 IF (azimuth > (90 - t) .AND. azimuth <= (270 - t))
THEN
881 IF (azimuth > (180 - t) .AND. azimuth <= (360 - t))
THEN
886 IF (azimuth <= (90 - t) .OR. azimuth > (270 - t))
THEN
893 svfbuveg = svfs - (1.0 - svfsveg)*(1.0 - psi)
894 gvfalb = shadow*alb_ground
895 gvfalbnosh = (1 - shadow)*alb_ground
896 kuptowall = (gvfalb*radi*sin(altitude*
deg2rad)) &
897 + (radd*svfbuveg + alb_bldg*(1 - svfbuveg)*(radg*(1 - f_sh) + radd*f_sh))*gvfalbnosh
899 CALL kvikt_veg(svfe, svfeveg, vikttot, viktveg, viktwall)
900 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
901 keastdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
902 keast = keasti + keastdg
904 CALL kvikt_veg(svfs, svfsveg, vikttot, viktveg, viktwall)
905 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
906 ksouthdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
907 ksouth = ksouthi + ksouthdg
909 CALL kvikt_veg(svfw, svfwveg, vikttot, viktveg, viktwall)
910 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
911 kwestdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
912 kwest = kwesti + kwestdg
914 CALL kvikt_veg(svfn, svfnveg, vikttot, viktveg, viktwall)
915 svfviktbuveg = (viktwall + (viktveg)*(1 - psi))
916 knorthdg = (radd*(1 - svfviktbuveg) + alb_bldg*(svfviktbuveg*(radg*(1 - f_sh) + radd*f_sh)) + kuptowall)*0.5
917 knorth = knorthi + knorthdg
926 REAL(KIND(1D0)),
INTENT(in) :: vikttot
927 REAL(KIND(1D0)),
INTENT(in) :: svf
928 REAL(KIND(1D0)),
INTENT(in) :: svfveg
929 REAL(KIND(1D0)),
INTENT(out) :: viktveg, viktwall
930 REAL(KIND(1D0)) :: svfvegbu
940 svfvegbu = (svfveg + svf - 1)
945 viktveg =
cal_vikt(svfvegbu, vikttot)
946 viktveg = viktveg - viktwall
951 REAL(kind(1d0)),
INTENT(in) :: svf_x, vikttot
952 REAL(kind(1d0)) :: vikt
955 - (63.227*svf_x**6 - 161.51*svf_x**5 &
956 + 156.91*svf_x**4 - 70.424*svf_x**3 &
957 + 16.773*svf_x**2 - 0.4863*svf_x))/vikttot
961 SUBROUTINE lwalls(svf, svfveg, svfaveg, &
963 altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, ldown, svfalfa, F_sh_in, & !input
964 Least, Lnorth, Lsouth, Lwest)
967 REAL(KIND(1D0)),
INTENT(in) :: altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, ldown
969 REAL(KIND(1D0)),
INTENT(in) :: svfalfa, svf, svfveg, svfaveg, F_sh_in
970 REAL(KIND(1D0)),
INTENT(in) :: Ldown2d, Lup2d
971 REAL(KIND(1D0)),
INTENT(out) :: Least, Lnorth, Lsouth, Lwest
973 REAL(KIND(1D0)) :: vikttot, aziE, aziN, aziS, aziW, c
975 REAL(KIND(1D0)) :: svfalfaE, svfalfaS, svfalfaW, svfalfaN
976 REAL(KIND(1D0)) :: alfaB, betaB, betasun
977 REAL(KIND(1D0)) :: Lground, Lrefl, Lsky, Lsky_allsky, Lveg, Lwallsh, Lwallsun
984 INTEGER,
PARAMETER :: SOLWEIG_ldown = 1
986 REAL(KIND(1D0)) :: viktveg, viktsky, viktrefl, viktwall
987 REAL(KIND(1D0)) :: svfE, svfS, svfW, svfN, svfEveg, svfSveg, svfWveg, svfNveg
988 REAL(KIND(1D0)) :: svfEaveg, svfSaveg, svfWaveg, svfNaveg, F_sh
1041 azin = azimuth - 90 + t
1042 azie = azimuth - 180 + t
1043 azis = azimuth - 270 + t
1045 f_sh = 2.*f_sh_in - 1.
1047 IF (solweig_ldown == 1)
THEN
1049 lsky_allsky = emis_sky*sbc*((ta + 273.15)**4)*(1 - c) + c*sbc*((ta + 273.15)**4)
1055 CALL lvikt_veg(svfe, svfeveg, svfeaveg, vikttot, &
1056 viktveg, viktsky, viktrefl, viktwall)
1058 IF (altitude > 0)
THEN
1059 alfab = atan(svfalfae)
1060 betab = atan(tan((svfalfae)*f_sh))
1061 betasun = ((alfab - betab)/2) + betab
1062 IF (azimuth > (180 - t) .AND. azimuth <= (360 - t))
THEN
1063 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azie*(
pi/180)))**4)* &
1064 viktwall*(1 - f_sh)*cos(betasun)*0.5
1065 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1068 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1072 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1074 lsky = ((svfe + svfeveg - 1)*lsky_allsky)*viktsky*0.5
1075 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1077 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1078 least = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1081 CALL lvikt_veg(svfs, svfsveg, svfsaveg, vikttot, &
1082 viktveg, viktsky, viktrefl, viktwall)
1084 IF (altitude > 0)
THEN
1085 alfab = atan(svfalfas)
1086 betab = atan(tan((svfalfas)*f_sh))
1087 betasun = ((alfab - betab)/2) + betab
1088 IF (azimuth <= (90 - t) .OR. azimuth > (270 - t))
THEN
1089 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azis*(
pi/180)))**4)* &
1090 viktwall*(1 - f_sh)*cos(betasun)*0.5
1091 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1094 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1098 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1100 lsky = ((svfs + svfsveg - 1)*lsky_allsky)*viktsky*0.5
1101 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1103 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1104 lsouth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1107 CALL lvikt_veg(svfw, svfwveg, svfwaveg, vikttot, &
1108 viktveg, viktsky, viktrefl, viktwall)
1110 IF (altitude > 0)
THEN
1111 alfab = atan(svfalfaw)
1112 betab = atan(tan((svfalfaw)*f_sh))
1113 betasun = ((alfab - betab)/2) + betab
1114 IF (azimuth > (360 - t) .OR. azimuth <= (180 - t))
THEN
1115 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(aziw*(
pi/180)))**4)* &
1116 viktwall*(1 - f_sh)*cos(betasun)*0.5
1117 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1120 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1124 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1126 lsky = ((svfw + svfwveg - 1)*lsky_allsky)*viktsky*0.5
1127 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1129 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1130 lwest = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1133 CALL lvikt_veg(svfn, svfnveg, svfnaveg, vikttot, &
1134 viktveg, viktsky, viktrefl, viktwall)
1136 IF (altitude > 0)
THEN
1137 alfab = atan(svfalfan)
1138 betab = atan(tan((svfalfan)*f_sh))
1139 betasun = ((alfab - betab)/2) + betab
1140 IF (azimuth > (90 - t) .AND. azimuth <= (270 - t))
THEN
1141 lwallsun = sbc*emis_wall*((ta + 273.15 + tw*sin(azin*(
pi/180)))**4)* &
1142 viktwall*(1 - f_sh)*cos(betasun)*0.5
1143 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*f_sh*0.5
1146 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1150 lwallsh = sbc*emis_wall*((ta + 273.15)**4)*viktwall*0.5
1152 lsky = ((svfn + svfnveg - 1)*lsky_allsky)*viktsky*0.5
1153 lveg = sbc*emis_wall*((ta + 273.15)**4)*viktveg*0.5
1155 lrefl = (ldown2d + lup2d)*(viktrefl)*(1 - emis_wall)*0.5
1156 lnorth = lsky + lwallsun + lwallsh + lveg + lground + lrefl
1189 SUBROUTINE lvikt_veg(isvf, isvfveg, isvfaveg, vikttot, & ! input
1190 viktveg, viktsky, viktrefl, viktwall)
1194 REAL(KIND(1D0)),
INTENT(in) :: vikttot
1195 REAL(KIND(1D0)),
INTENT(in) :: isvf
1196 REAL(KIND(1D0)),
INTENT(in) :: isvfveg
1197 REAL(KIND(1D0)),
INTENT(in) :: isvfaveg
1198 REAL(KIND(1D0)),
INTENT(out) :: viktveg
1199 REAL(KIND(1D0)),
INTENT(out) :: viktsky
1200 REAL(KIND(1D0)),
INTENT(out) :: viktrefl
1201 REAL(KIND(1D0)),
INTENT(out) :: viktwall
1203 REAL(KIND(1D0)) :: viktonlywall
1204 REAL(KIND(1D0)) :: viktaveg
1205 REAL(KIND(1D0)) :: svfvegbu
1207 viktonlywall = (vikttot - &
1208 (63.227*isvf**6 - 161.51*isvf**5 + 156.91*isvf**4 &
1209 - 70.424*isvf**3 + 16.773*isvf**2 - 0.4863*isvf))/vikttot
1211 viktaveg = (vikttot &
1212 - (63.227*isvfaveg**6 - 161.51*isvfaveg**5 &
1213 + 156.91*isvfaveg**4 - 70.424*isvfaveg**3 &
1214 + 16.773*isvfaveg**2 - 0.4863*isvfaveg))/vikttot
1216 viktwall = viktonlywall - viktaveg
1218 svfvegbu = (isvfveg + isvf - 1)
1219 viktsky = (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1220 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1221 + 16.773*svfvegbu**2 - 0.4863*svfvegbu)/vikttot
1222 viktrefl = (vikttot &
1223 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1224 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1225 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1226 viktveg = (vikttot &
1227 - (63.227*svfvegbu**6 - 161.51*svfvegbu**5 &
1228 + 156.91*svfvegbu**4 - 70.424*svfvegbu**3 &
1229 + 16.773*svfvegbu**2 - 0.4863*svfvegbu))/vikttot
1230 viktveg = viktveg - viktwall
1237 REAL(KIND(1D0)) IX, MAXPOS, ISIGNM
1239 IF (ix < 0 .OR. ix > maxpos) isignm = -1
1240 IF (ix == 0) isignm = 0
1268 INTEGER,
INTENT(in) :: b
1269 INTEGER,
INTENT(out) :: mb
1270 INTEGER,
INTENT(out) :: md
1271 INTEGER,
INTENT(out) :: seas
1272 INTEGER,
INTENT(in) :: year
1273 INTEGER :: t1, t2, t3
1276 REAL(KIND(1D0)) :: latitude
1287 IF ((modulo(year, t1) == 0) .AND. (modulo(year, t2) /= 0) .OR. (modulo(year, t3) == 0))
THEN
1296 ELSEIF (b > 31 .AND. b <= 59 + k)
THEN
1299 ELSEIF (b > 59 + k .AND. b <= 90 + k)
THEN
1302 ELSEIF (b > 90 + k .AND. b <= 120 + k)
THEN
1305 ELSEIF (b > 120 + k .AND. b <= 151 + k)
THEN
1308 ELSEIF (b > 151 + k .AND. b <= 181 + k)
THEN
1311 ELSEIF (b > 181 + k .AND. b <= 212 + k)
THEN
1314 ELSEIF (b > 212 + k .AND. b <= 243 + k)
THEN
1317 ELSEIF (b > 243 + k .AND. b <= 273 + k)
THEN
1320 ELSEIF (b > 273 + k .AND. b <= 304 + k)
THEN
1323 ELSEIF (b > 304 + k .AND. b <= 334 + k)
THEN
1326 ELSEIF (b > 334 + k)
THEN
1332 IF (latitude > 0)
THEN
1333 IF (mb > 3 .AND. mb < 10)
THEN
1339 IF (mb < 4 .OR. mb > 9)
THEN
1350 INTEGER :: mon, ne, k, b
1354 ELSE IF (mon == 2)
THEN
1356 ELSE IF (mon == 3)
THEN
1358 ELSE IF (mon == 4)
THEN
1360 ELSE IF (mon == 5)
THEN
1362 ELSE IF (mon == 6)
THEN
1364 ELSE IF (mon == 7)
THEN
1366 ELSE IF (mon == 8)
THEN
1369 ELSE IF (mon == 9)
THEN
1371 ELSE IF (mon == 10)
THEN
1373 ELSE IF (mon == 11)
THEN
1375 ELSE IF (mon == 12)
THEN
1385 INTEGER :: nroDays, year_int
1387 IF (mod(year_int, 100) /= 0 .AND. mod(year_int, 4) == 0)
THEN
1389 ELSEIF (mod(year_int, 400) == 0)
THEN
1401 INTEGER,
INTENT(in) :: year_int
1404 IF (mod(year_int, 100) /= 0 .AND. mod(year_int, 4) == 0)
THEN
1406 ELSEIF (mod(year_int, 400) == 0)
THEN
1423 INTEGER DATE, MONTH, DAY, YR, MN, N1, N2, DOW, YEAR
1431 IF (mn > 2)
GO TO 10
143410 n1 = (26*(mn + 1))/10
1436 day = (date + n1 + n2 - (yr/100) + (yr/400) - 1)
1437 dow = mod(day, 7) + 1
1448 INTEGER :: HOURS, MINS, doy
1449 REAL(KIND(1D0)) :: dectime, SECS, DH, DM, DS
1452 doy = floor(dectime)
1469 SUBROUTINE daylen(DOY, XLAT, DAYL, DEC, SNDN, SNUP)
1477 REAL(KIND(1D0)),
INTENT(IN) :: XLAT
1478 REAL(KIND(1D0)),
INTENT(OUT) :: DEC, DAYL, SNDN, SNUP
1479 REAL(KIND(1D0)) :: SOC
1480 REAL(KIND(1D0)),
PARAMETER :: RAD =
pi/180.0
1485 dec = -23.45*cos(2.0*
pi*(doy + 10.0)/365.0)
1488 soc = tan(rad*dec)*tan(rad*xlat)
1489 soc = min(max(soc, -1.0), 1.0)
1492 dayl = 12.0 + 24.0*asin(soc)/
pi
1493 snup = 12.0 - dayl/2.0
1494 sndn = 12.0 + dayl/2.0
1514 id, it, imin, isec, & ! input
1517 INTEGER,
INTENT(in) :: id, it, imin, isec
1519 REAL(KIND(1D0)),
INTENT(out) :: dectime
1521 dectime = real(id - 1, kind(1d0)) &
1522 + real(it, kind(1d0))/24 &
1523 + real(imin, kind(1d0))/(60*24) &
1524 + real(isec, kind(1d0))/(60*60*24)
1531 nsh, nsh_real, tstep_real)
1533 INTEGER,
INTENT(in) :: tstep
1535 INTEGER,
INTENT(out) :: nsh
1536 REAL(KIND(1D0)),
INTENT(out) :: nsh_real
1537 REAL(KIND(1D0)),
INTENT(out) :: tstep_real
1540 tstep_real = tstep*1.0
1545 iy, id, lat, & !input
1549 INTEGER,
INTENT(in) :: iy
1550 INTEGER,
INTENT(in) :: id
1551 REAL(KIND(1D0)),
INTENT(in) :: lat
1553 INTEGER,
DIMENSION(3),
INTENT(OUT) :: dayofWeek_id
1560 CALL day2month(id, mb, date, seas, iy, lat)
1563 dayofweek_id(1) = wd
1564 dayofweek_id(2) = mb
1565 dayofweek_id(3) = seas
1570 id, startDLS, endDLS, & !input
1574 INTEGER,
INTENT(in) :: id, startDLS, endDLS
1575 INTEGER,
INTENT(out) :: DLS
1578 IF (id > startdls .AND. id < enddls) dls = 1
real(kind(1d0)), parameter deg2rad
integer, parameter ncolumnsdataoutbeers
real(kind(1d0)), parameter rad2deg
subroutine suews_cal_tstep(tstep, nsh, nsh_real, tstep_real)
real(kind(1d0)) function cal_vikt(svf_x, vikttot)
subroutine sun_distance(jday, D)
subroutine tsurfbeers(iy, Ta, RH, radI, I0, dectime, SNUP, altitude, zen, timezone, lat, lng, alt, Tg, Tgwall, altmax)
subroutine shadowgroundkusaka(HW, azimuth, zen, shadowground, shadowwalls)
subroutine leapyearcalc(year_int, nroDays)
subroutine cal_ci_latenight(iy, DOY, Ta_degC, RH_frac, radG, lat, P_kPa, CIlatenight, dectime_sunrise, zen_sunrise, I0_sunrise)
subroutine suews_cal_weekday(iy, id, lat, dayofWeek_id)
subroutine daylen(DOY, XLAT, DAYL, DEC, SNDN, SNUP)
subroutine kroof(radI, radD, radG, F_sh, altitude, svfr, svfveg, shadow, psi, alb_bldg, Kdown)
subroutine month2day(mon, ne, k, b)
subroutine cal_svfalfa(svfr, svfveg, svfalfa, tmp)
subroutine cylindric_wedge(zen, svfalfa, F_sh)
subroutine lvikt_veg(isvf, isvfveg, isvfaveg, vikttot, viktveg, viktsky, viktrefl, viktwall)
subroutine day_of_week(DATE, MONTH, YEAR, DOW)
real(kind(1d0)), parameter pi
subroutine dectime_to_timevec(dectime, HOURS, MINS, SECS)
subroutine kwalls(svf, svfveg, shadow, F_sh, radI, radG, radD, azimuth, altitude, psi, t, alb_ground, alb_bldg, Keast, Knorth, Ksouth, Kwest)
subroutine suews_cal_dls(id, startDLS, endDLS, DLS)
real(kind(1d0)) function hwtosvf_ground(hw)
subroutine beers_cal_main(iy, id, dectime, lamdaP, lamdaF, avkdn, ldown, Temp_C, avrh, Press_hPa, Tsurf, lat, lng, alt, timezone, zenith_deg, azimuth, alb_ground, alb_bldg, emis_ground, emis_wall, dataOutLineBEERS)
elemental integer function days_of_year(year_int)
subroutine issign(IX, MAXPOS, ISIGNM)
subroutine lwalls(svf, svfveg, svfaveg, Ldown2d, Lup2d, altitude, Ta, Tw, SBC, emis_wall, emis_sky, t, CI, azimuth, ldown, svfalfa, F_sh_in, Least, Lnorth, Lsouth, Lwest)
real(kind(1d0)) function hwtosvf_roof(hw)
subroutine diffusefraction(radG, altitude, Kt, Ta, RH, radI, radD)
subroutine suews_cal_dectime(id, it, imin, isec, dectime)
subroutine day2month(b, mb, md, seas, year, latitude)
real(kind(1d0)) function cal_ratio_height2width(lamdaP, lamdaF)
subroutine clearnessindex_2013b(zen, DOY, Ta_degC, RH_frac, radG, lat, P_kPa, I0, CI, Kt, I0et, CIuncorr)
subroutine kvikt_veg(svf, svfveg, vikttot, viktveg, viktwall)
subroutine narp_cal_sunposition(year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)