SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
solweig_module Module Reference

Functions/Subroutines

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)
 
real(kind(1d0)) function cal_ratio_height2width (lamdaP, lamdaF)
 
real(kind(1d0)) function hwtosvf_ground (hw)
 
real(kind(1d0)) function hwtosvf_roof (hw)
 
subroutine clearnessindex_2013b (zen, DOY, Ta, RH, radG, lat, P_kPa, I0, CI, Kt, I0et, CIuncorr)
 
subroutine sun_distance (jday, D)
 
subroutine cylindric_wedge (zen, svfalfa, F_sh)
 
subroutine diffusefraction (radG, altitude, Kt, Ta, RH, radI, radD)
 
subroutine kside_veg_v24 (shadow, F_sh, radI, radG, radD, azimuth, altitude, psi, t, albedo, Keast, Knorth, Ksouth, Kwest)
 
subroutine kvikt_veg (isvf, isvfveg, vikttot, viktveg, viktwall)
 
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 lvikt_veg (isvf, isvfveg, isvfaveg, vikttot, viktveg, viktsky, viktrefl, viktwall)
 
subroutine issign (IX, MAXPOS, ISIGNM)
 
character(len=20) function str (k)
 
subroutine shadowingfunction_urban (azimuth, altitude, scale, shadow)
 
subroutine sunonsurface_veg (iazimuthA, scale, buildings, first, second, psi, sos)
 

Function/Subroutine Documentation

◆ cal_ratio_height2width()

real(kind(1d0)) function solweig_module::cal_ratio_height2width ( real(kind(1d0)), intent(in)  lamdaP,
real(kind(1d0)), intent(in)  lamdaF 
)

Definition at line 460 of file suews_phys_solweig.f95.

Referenced by solweig_cal_main().

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 
real(kind(1d0)) hw
Here is the caller graph for this function:

◆ clearnessindex_2013b()

subroutine solweig_module::clearnessindex_2013b ( real(kind(1d0)), intent(in)  zen,
integer, intent(in)  DOY,
real(kind(1d0)), intent(in)  Ta,
real(kind(1d0)), intent(in)  RH,
real(kind(1d0)), intent(in)  radG,
real(kind(1d0)), intent(in)  lat,
real(kind(1d0)), intent(in)  P_kPa,
real(kind(1d0)), intent(out)  I0,
real(kind(1d0)), intent(out)  CI,
real(kind(1d0)), intent(out)  Kt,
real(kind(1d0)), intent(out)  I0et,
real(kind(1d0)), intent(out)  CIuncorr 
)

Definition at line 508 of file suews_phys_solweig.f95.

References sun_distance().

Referenced by solweig_cal_main().

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
real(kind(1d0)) a2
real(kind(1d0)), parameter pi
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cylindric_wedge()

subroutine solweig_module::cylindric_wedge ( real(kind(1d0)), intent(in)  zen,
real(kind(1d0)), dimension(1, 1), intent(in)  svfalfa,
real(kind(1d0)), dimension(1, 1), intent(out)  F_sh 
)

Definition at line 638 of file suews_phys_solweig.f95.

Referenced by lside_veg_v2(), and solweig_cal_main().

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 
real(kind(1d0)), parameter pi
Here is the caller graph for this function:

◆ diffusefraction()

subroutine solweig_module::diffusefraction ( real(kind(1d0)), intent(in)  radG,
real(kind(1d0)), intent(in)  altitude,
real(kind(1d0)), intent(in)  Kt,
real(kind(1d0)), intent(in)  Ta,
real(kind(1d0)), intent(in)  RH,
real(kind(1d0)), intent(out)  radI,
real(kind(1d0)), intent(out)  radD 
)

Definition at line 721 of file suews_phys_solweig.f95.

References allocatearray::deg2rad.

Referenced by solweig_cal_main().

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
real(kind(1d0)), parameter deg2rad
Here is the caller graph for this function:

◆ hwtosvf_ground()

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

Definition at line 478 of file suews_phys_solweig.f95.

Referenced by solweig_cal_main().

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 
real(kind(1d0)) hw
Here is the caller graph for this function:

◆ hwtosvf_roof()

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

Definition at line 493 of file suews_phys_solweig.f95.

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 
real(kind(1d0)) hw

◆ issign()

subroutine solweig_module::issign ( real(kind(1d0))  IX,
real(kind(1d0))  MAXPOS,
real(kind(1d0))  ISIGNM 
)

Definition at line 1425 of file suews_phys_solweig.f95.

Referenced by shadowingfunction_urban(), and sunonsurface_veg().

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
Here is the caller graph for this function:

◆ kside_veg_v24()

subroutine solweig_module::kside_veg_v24 ( real(kind(1d0)), dimension(1, 1), intent(in)  shadow,
real(kind(1d0)), dimension(1, 1), intent(in)  F_sh,
real(kind(1d0)), intent(in)  radI,
real(kind(1d0)), intent(in)  radG,
real(kind(1d0)), intent(in)  radD,
real(kind(1d0)), intent(in)  azimuth,
real(kind(1d0)), intent(in)  altitude,
real(kind(1d0)), intent(in)  psi,
real(kind(1d0)), intent(in)  t,
real(kind(1d0)), intent(in)  albedo,
real(kind(1d0)), dimension(1, 1), intent(out)  Keast,
real(kind(1d0)), dimension(1, 1), intent(out)  Knorth,
real(kind(1d0)), dimension(1, 1), intent(out)  Ksouth,
real(kind(1d0)), dimension(1, 1), intent(out)  Kwest 
)

Definition at line 1049 of file suews_phys_solweig.f95.

References kvikt_veg().

Referenced by solweig_cal_main().

1049 
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)
real(kind(1d0)), parameter pi
Here is the call graph for this function:
Here is the caller graph for this function:

◆ kvikt_veg()

subroutine solweig_module::kvikt_veg ( real(kind(1d0)), dimension(1, 1), intent(in)  isvf,
real(kind(1d0)), dimension(1, 1), intent(in)  isvfveg,
real(kind(1d0)), intent(in)  vikttot,
real(kind(1d0)), dimension(1, 1), intent(out)  viktveg,
real(kind(1d0)), dimension(1, 1), intent(out)  viktwall 
)

Definition at line 1131 of file suews_phys_solweig.f95.

Referenced by kside_veg_v24().

1131 
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
Here is the caller graph for this function:

◆ lside_veg_v2()

subroutine solweig_module::lside_veg_v2 ( real(kind(1d0)), dimension(1, 1), intent(in)  Ldown2d,
real(kind(1d0)), dimension(1, 1), intent(in)  Lup2d,
real(kind(1d0)), intent(in)  altitude,
real(kind(1d0)), intent(in)  Ta,
real(kind(1d0)), intent(in)  Tw,
real(kind(1d0)), intent(in)  SBC,
real(kind(1d0)), intent(in)  emis_wall,
real(kind(1d0)), intent(in)  emis_sky,
real(kind(1d0)), intent(in)  t,
real(kind(1d0)), intent(in)  CI,
real(kind(1d0)), intent(in)  azimuth,
real(kind(1d0)), intent(in)  zen,
real(kind(1d0)), intent(in)  ldown,
real(kind(1d0)), dimension(1, 1), intent(in)  svfalfa,
real(kind(1d0)), dimension(1, 1), intent(out)  Least,
real(kind(1d0)), dimension(1, 1), intent(out)  Lnorth,
real(kind(1d0)), dimension(1, 1), intent(out)  Lsouth,
real(kind(1d0)), dimension(1, 1), intent(out)  Lwest 
)

Definition at line 1157 of file suews_phys_solweig.f95.

References cylindric_wedge(), and lvikt_veg().

Referenced by solweig_cal_main().

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 
real(kind(1d0)), parameter pi
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lvikt_veg()

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

Definition at line 1373 of file suews_phys_solweig.f95.

Referenced by lside_veg_v2().

1373 
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 
Here is the caller graph for this function:

◆ shadowingfunction_urban()

subroutine solweig_module::shadowingfunction_urban ( real(kind(1d0)), intent(in)  azimuth,
real(kind(1d0)), intent(in)  altitude,
real(kind(1d0)), intent(in)  scale,
real(kind(1d0)), dimension(1, 1), intent(out)  shadow 
)

Definition at line 1523 of file suews_phys_solweig.f95.

References issign().

Referenced by solweig_cal_main().

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 
real(kind(1d0)), parameter pi
Here is the call graph for this function:
Here is the caller graph for this function:

◆ solweig_cal_main()

subroutine solweig_module::solweig_cal_main ( integer, intent(in)  id,
integer, intent(in)  it,
real(kind(1d0)), intent(in)  dectime,
real(kind(1d0)), intent(in)  lamdaP,
real(kind(1d0)), intent(in)  lamdaF,
real(kind(1d0)), intent(in)  avkdn,
real(kind(1d0)), intent(in)  ldown,
real(kind(1d0)), intent(in)  Temp_C,
real(kind(1d0)), intent(in)  avrh,
real(kind(1d0)), intent(in)  Press_hPa,
real(kind(1d0)), intent(in)  Tg,
real(kind(1d0)), intent(in)  lat,
real(kind(1d0)), intent(in)  zenith_deg,
real(kind(1d0)), intent(in)  azimuth,
real(kind(1d0)), intent(in)  scale,
real(kind(1d0)), intent(in)  alb_ground,
real(kind(1d0)), intent(in)  alb_bldg,
real(kind(1d0)), intent(in)  emis_ground,
real(kind(1d0)), intent(in)  emis_wall,
real(kind(1d0)), intent(in)  heightgravity,
real(kind(1d0)), dimension(ncolumnsdataoutsol - 5), intent(out)  dataOutLineSOLWEIG 
)

Definition at line 75 of file suews_phys_solweig.f95.

References cal_ratio_height2width(), clearnessindex_2013b(), cylindric_wedge(), daylen(), allocatearray::deg2rad, diffusefraction(), hwtosvf_ground(), kside_veg_v24(), lside_veg_v2(), shadowingfunction_urban(), and sunonsurface_veg().

Referenced by suews_driver::suews_cal_main().

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 
real(kind(1d0)) hw
real(kind(1d0)), parameter deg2rad
integer id
real(kind(1d0)) alb_ground
integer it
real(kind(1d0)) dectime
real(kind(1d0)), parameter pi
subroutine daylen(DOY, XLAT, DAYL, DEC, SNDN, SNUP)
real(kind(1d0)), dimension(ncolumnsdataoutsol - 5) dataoutlinesolweig
Here is the call graph for this function:
Here is the caller graph for this function:

◆ str()

character(len=20) function solweig_module::str ( integer, intent(in)  k)

Definition at line 1435 of file suews_phys_solweig.f95.

1435  INTEGER, INTENT(in) :: k
1436  WRITE (str, *) k
1437  str = adjustl(str)
real(kind(1d0)) k

◆ sun_distance()

subroutine solweig_module::sun_distance ( integer  jday,
real(kind(1d0))  D 
)

Definition at line 624 of file suews_phys_solweig.f95.

Referenced by clearnessindex_2013b().

624 
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 
Here is the caller graph for this function:

◆ sunonsurface_veg()

subroutine solweig_module::sunonsurface_veg ( real(kind(1d0)), intent(in)  iazimuthA,
real(kind(1d0)), intent(in)  scale,
real(kind(1d0)), dimension(1, 1), intent(in)  buildings,
integer  first,
integer  second,
real(kind(1d0)), intent(in)  psi,
real(kind(1d0)), dimension(1, 1), intent(out)  sos 
)

Definition at line 1861 of file suews_phys_solweig.f95.

References issign().

Referenced by solweig_cal_main().

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 
real(kind(1d0)), parameter pi
Here is the call graph for this function:
Here is the caller graph for this function: