!! Manage solar radiation with arbitrary time cumulation 
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.1 - 16th June 2021  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 19/May/2021 | Original code |
! | 1.1      | 26/Jun/2021 | added net radiation computation |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Module to manage solar radiation with arbitrary time cumulation
!
MODULE SolarRadiation


! Modules used: 

USE DataTypeSizes, ONLY : &  
  ! Imported Parameters:
  float, short


USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE ObservationalNetworks, ONLY : &
  !Imported types:
  ObservationalNetwork, &
  !Imported routines:
  ReadMetadata, ReadData, &
  CopyNetwork, &
  !Imported operators:
  ASSIGNMENT (=)

USE Chronos, ONLY : &
  !Imported types:
  DateTime, &
  !Imported variables:
  timeString, &
  !imported routines:
  DayOfYear, GetHour, &
  !Imported operations:
  OPERATOR (<), &
  OPERATOR (>), &
  OPERATOR (>=), &
  OPERATOR (<=), &
  OPERATOR (+), &
  OPERATOR (==), &
  ASSIGNMENT (=)

USE GridLib, ONLY : &
  !Imported types:
  grid_real, grid_integer, &
  !Imported routines:
  NewGrid, GetDtGrid, &
  GridDestroy, ExportGrid, &
  SetLongName, SetStandardName, &
  SetUnits, SetCurrentTime, &
  !Imported parameters:
  NET_CDF, ESRI_ASCII, &
  ESRI_BINARY
  

USE GridOperations, ONLY : &
  !Imported routines:
  GridByIni, GridResample, &
  CRSisEqual, GridConvert, &
  GetXY, GetIJ, &
  !Imported operations:
  ASSIGNMENT (=)

USE Units, ONLY : &
    !imported parameters:
    degToRad, pi, hour, stefanBoltzman

USE IniLib, ONLY : &
  !Imported types:
  IniList, &
  !Imported routines:
  IniReadInt, IniReadReal, &
  KeyIsPresent, IniReadString

USE Utilities, ONLY: &
  !Imported routines:
  GetUnit

USE MeteoUtilities, ONLY: &
  !Imported routines:
  ReadField, Interpolate, &
  Offset, Scale

USE FileSys, ONLY: &
  !Imported routines:
  DirExists, FileExists, &
  FileDelete, CurrentDir, &
  GetOS, DirNew, &
  !Imported parameters:
  WIN32

USE Morphology, ONLY: &
    !imported routines:
    DeriveAspect, DeriveSlope
    
USE GeoLib, ONLY: &
  !imported type
  Coordinate, &
  !Imported variables:
  point1, point2, &
  !Imported routines:
  DecodeEPSG, Convert, &
  Distance, &
  !Imported types:
  CRS, &
  !Imported operations:
  OPERATOR (==)

USE DomainProperties, ONLY : &
    !imported variables
    latCentroid 


USE StringManipulation, ONLY: &
   !Imported routines:
   ReplaceChar, ToString
    
IMPLICIT NONE 

! Global (i.e. public) Declarations: 

INTEGER (KIND = short)     :: dtRadiation  = 0!!cumulation deltat 
TYPE (ObservationalNetwork) :: radiometers !!radiation stations network
TYPE (grid_real)            :: radiation	   	  !![W/m2]
TYPE (grid_real)            :: netRadiation !!incoming and outcoming shortwave and longwave radiation [W/m2]
TYPE (grid_real)            :: netRadiationFAO !!net radiation of FAO reference grass with albedo = 0.23


!Private declarations

!define type for viewing angle: the maximum angle of terrain obstruction
TYPE ViewingAngle
     REAL (KIND = float) :: angle (16)
END TYPE ViewingAngle



TYPE(DateTime), PRIVATE         :: timeNew !!time when new  data must be read
INTEGER (KIND = short), PRIVATE :: fileunit !! unit of file containing  data
INTEGER (KIND = short), PRIVATE :: interpolationMethod  !!method to spatial interpolate site data
INTEGER (KIND = short), PRIVATE :: interpolationMethod_assignment  !!method to assign spatial interpolation !
                                                                   !!1 = one method for the entire domain, 2 = a map with interpolation method codes
TYPE (grid_integer),    PRIVATE :: interpolationMethod_map
INTEGER (KIND = short), PRIVATE :: interpolationMethod_vector (3) !!defines active interpolation methods 
TYPE (grid_real),       PRIVATE :: interpolatedMap (3)  !!1 = map for thiessen, 2 = map for idw, 3 = map for kriging
REAL (KIND = float),    PRIVATE :: cellsizeInterpolation !!spatial resolution of interpolated grid
INTEGER (KIND = short), PRIVATE :: neighbors  !!number of closest data to use for interpolation
TYPE(grid_real),        PRIVATE :: grid_devst !!standard deviation of kriging interpolation
CHARACTER (LEN = 300),  PRIVATE :: fileGrid !!file containing grids
INTEGER (KIND = short), PRIVATE :: dtGrid !!dt of imported  field
INTEGER (KIND = short), PRIVATE :: elevationDrift !! 1 = use elevation to modify interpolated data
INTEGER (KIND = short), PRIVATE :: timezone !!local timezone
TYPE(grid_real),        PRIVATE :: shadowGrid !!shadow grid 1 = shadow, 0 = no shadow
TYPE (ViewingAngle), ALLOCATABLE, PRIVATE :: viewangle (:,:) !!raster of maximum viewing angles
REAL (KIND = float),   PRIVATE :: azimuth (16) 
TYPE(grid_real),        PRIVATE :: slope !!terreain slope to be used to modify interpolated data [radians]
TYPE(grid_real),        PRIVATE :: aspect !!terreain aspect to be used to modify interpolated data [radians]
INTEGER (KIND = short), PRIVATE :: i, j  !!loop index
REAL (KIND = float), PRIVATE, PARAMETER :: MISSING_DEF_REAL = -9999.9
REAL (KIND = float), PRIVATE  :: scale_factor  !!scale factor to apply to final map
REAL (KIND = float), PRIVATE  :: offset_value !!offset to apply to final map
REAL (KIND = float), PRIVATE :: valid_prcn !!when data from several time steps are read, 
                                      !!this is the minimum percentage (0-1) of valid 
                                       !!data that must be prresent to consider 
                                      !!valid the aggregated value.

!used by inverse distance weighting interpolation
REAL (KIND = float), PRIVATE :: idw_power !! power to be used with IDW

!used by kriging
INTEGER (KIND = short), PRIVATE :: krige_var !! when set to 1 a grid of kriging variance is generated and exported if export option is active, default to 0
INTEGER (KIND = short), PRIVATE :: krige_anisotropy !! considers anisotropy, default = 0 excludes anisotropy
INTEGER (KIND = short), PRIVATE :: krige_varmodel !! 1 = spherical, 2 = exponential, 3 = gaussian, 0 = automatic fitting. default to 2
INTEGER (KIND = short), PRIVATE :: krige_lags !! number of lags for semivariogram. if undefined or set to 0 default to 15
REAL (KIND = float)   , PRIVATE :: krige_maxlag !! maximum distance to be considered for semivariogram assessment. If undefined or set to 0, it is computed automatically

INTEGER (KIND = short), PRIVATE :: export !!activates grid exporting
CHARACTER (LEN = 1000), PRIVATE :: export_path  !!folder where to put exported grids
CHARACTER (LEN = 1000), PRIVATE :: export_file  !!name of  exported  file
CHARACTER (LEN = 1000), PRIVATE :: export_file_net  !!name of  exported  net radiation file
CHARACTER (LEN = 1000), PRIVATE :: export_file_var  !!name of  exported  variance grid
TYPE (DateTime), PRIVATE :: export_start   !! time and date to start exporting
TYPE (DateTime), PRIVATE :: export_stop !! time and date to stop exporting
INTEGER (KIND = short), PRIVATE ::  export_dt !! time between two exportations
TYPE (DateTime), PRIVATE :: timeNewExport !! time when new exporting must occur
INTEGER (KIND = short), PRIVATE :: export_format !! 1 = esri_ascii, 2 = esri_binary, 3 = netcdf 
INTEGER (KIND = short), PRIVATE :: export_epsg !! coordinate reference system of exported grid 
TYPE (grid_real), PRIVATE :: exportedGrid, exportedGridVar, exportedGridNet
TYPE (CRS), PRIVATE :: exportGridMapping
LOGICAL, PRIVATE   :: needConversion



!Public routines
PUBLIC :: SolarRadiationInit
PUBLIC :: SolarRadiationRead


!Local routines
PRIVATE :: SetSpecificProperties
PRIVATE :: SkyView
PRIVATE :: CastShadow
PRIVATE :: SolarDeclination
PRIVATE :: SunElevationAngle
PRIVATE :: SolarHourAngle
PRIVATE :: SolarAzimuth
PRIVATE :: OpticalDepth
PRIVATE :: ComputeNetRadiation


!=======
CONTAINS
!=======
! Define procedures contained in this module.

!==============================================================================
!| Description:
!   Initialize solar radiation
SUBROUTINE SolarRadiationInit &
!
( ini, mask, dtMeteo, tstart, dem, dem_loaded, albedo_loaded, &
  dtTemperature, dtRelHumidity )

IMPLICIT NONE

TYPE (IniList), INTENT(IN) :: ini
TYPE (grid_integer), INTENT(IN) :: mask  !!defines interpolation extent
INTEGER (KIND = short), INTENT(IN) :: dtMeteo !! deltat of meteo data reading
TYPE (DateTime),     INTENT(IN) :: tstart !!initial time
TYPE(grid_real),    INTENT(in)  :: dem !!digital elevation model to be used to modify interpolated data
LOGICAL , INTENT (in) :: dem_loaded !! true if dem has been loaded
LOGICAL , INTENT (in) :: albedo_loaded !! true if dem has been loaded
INTEGER (KIND = short), INTENT(IN) :: dtTemperature !! delta time temperature
INTEGER (KIND = short), INTENT(IN) :: dtRelHumidity !! delta time relative humidity



!local declarations
CHARACTER (LEN = 1000) :: filename
TYPE (grid_real) :: meteoTemp


!-------------------------end of declarations----------------------------------
  
!check if albedo and dem have been loaded by domain properties
IF ( .NOT. dem_loaded) THEN
    CALL Catch ('error', 'SolarRadiation',   &
                'dem  was not loaded ')
END IF
      
IF ( .NOT. albedo_loaded) THEN
    CALL Catch ('error', 'SolarRadiation',   &
                'albedo  was not loaded ')
END IF

!check if air temperature and relative humidity have been initialized
IF ( dtTemperature <= 0) THEN
    CALL Catch ('error', 'SolarRadiation',   &
                'air temperature not initialized ')
END IF
      
IF ( dtRelHumidity <= 0) THEN
    CALL Catch ('error', 'SolarRadiation',   &
                'air relative humidity not initialized ')
END IF

  !set initial time
  timeNew = tstart  

  !initialize grid
  CALL NewGrid (radiation, mask, 0.)
  CALL NewGrid (grid_devst, mask, 0.)
  
  !set deltat
  dtRadiation = IniReadInt ('dt', ini, section = 'solar-radiation')

	!check dt is multiple of dtMeteo
  IF (.NOT.(MOD(dtRadiation,dtMeteo) == 0)) THEN
           CALL Catch ('error', 'SolarRadiation',   &
		                 'dt is not multiple of dtMeteo')
  END IF
  
  !set valid threshold
  IF (KeyIsPresent ('valid-threshold', ini, &
                         section = 'solar-radiation') ) THEN
      valid_prcn = IniReadReal ('valid-threshold', ini, &
                                        section = 'solar-radiation')
  ELSE ! set to default = 1.0 
      CALL Catch ('info', 'SolarRadiation',   &
		                 'valid-threshold not defined, set to 1')
      valid_prcn = 1.
  END IF

  !set cell-size
  cellsizeInterpolation = mask % cellsize

  
  !interpolation-assignment method
  IF (KeyIsPresent('interpolation-assignment', ini, &
                     section = 'solar-radiation') ) THEN
      interpolationMethod_assignment = IniReadInt ('interpolation-assignment', &
                                     ini, section = 'solar-radiation')
  ELSE
    CALL Catch ('error', 'SolarRadiation',   &
                'interpolation-assignment missing in meteo configuration file')
  END IF   

  !set interpolation method
  IF (interpolationMethod_assignment == 1) THEN

      interpolationMethod = IniReadInt ('interpolation', ini, &
                                    section = 'solar-radiation')
      
	    CALL SetSpecificProperties (interpolationMethod, ini, mask)
  
  ELSE !read map
      interpolationMethod = -1
      CALL GridByIni (ini, interpolationMethod_map, 'solar-radiation', &
                     "interpolation")
      !scan for interpolation methods included
      interpolationMethod_vector = 0
      DO i = 1, interpolationMethod_map % idim
        DO j = 1, interpolationMethod_map % jdim
            IF (interpolationMethod_map % mat(i,j) /= &
                interpolationMethod_map % nodata) THEN
                interpolationMethod_vector &
                   (interpolationMethod_map % mat (i,j)) = &
                    interpolationMethod_map % mat (i,j)
            END IF
        END DO
      END DO

      DO i = 1, 3
         CALL SetSpecificProperties (interpolationMethod_vector (i), ini, mask)
      END DO

  END IF
  
  
  !allocate variable for netradiation
  CALL NewGrid (netRadiation, mask) 
  CALL NewGrid (netRadiationFAO, mask) 
  
  !scale factor and offset
  IF (KeyIsPresent('offset', ini, section = 'solar-radiation')) THEN	
	offset_value = IniReadReal ('offset', ini, section = 'solar-radiation')
  ELSE
	offset_value = MISSING_DEF_REAL
  END IF
	
  IF (KeyIsPresent('scale-factor', ini, section = 'solar-radiation')) THEN	
	scale_factor = IniReadReal ('scale-factor', ini, section = 'solar-radiation')
  ELSE
	scale_factor = MISSING_DEF_REAL
  END IF			

   !set power_idw
   IF (KeyIsPresent('idw-power', ini, section = 'solar-radiation')) THEN	
	   idw_power = IniReadReal ('idw-power', ini, section = 'solar-radiation')
   ELSE !set default value
	    idw_power = 2.
   END IF	
 
	!file
    filename = IniReadString ('file', ini, section = 'solar-radiation')
    IF (interpolationMethod_assignment == 1 .AND. &
        interpolationMethod == 0) THEN !data are stored in net-cdf file
        
       !store net-cdf filename
       fileGrid = filename
       
       IF ( KeyIsPresent ('variable', ini, section = 'solar-radiation') ) THEN
          radiation % var_name = IniReadString ('variable', &
                                             ini, section = 'solar-radiation')
            !read grid and store as temporary file
             CALL NewGrid (layer = meteoTemp, fileName = TRIM(fileGrid), &
                         fileFormat = NET_CDF, &
                         variable = TRIM(radiation % var_name) )
       ELSE IF  (KeyIsPresent ('standard_name', ini, &
                              section = 'solar-radiation') ) THEN
          radiation % standard_name = IniReadString ('standard_name', &
                                              ini, section = 'solar-radiation')
       ELSE
          CALL Catch ('error', 'SolarRadiation',   &
		       'standard_name or variable missing in section temeprature' )
       END IF      

       !set cellsize to zero. Optimal cellsize is automatically computed
       cellsizeInterpolation = 0.
            
    ELSE !open file containing local measurements
       fileunit = GetUnit ()
	     OPEN(fileunit, file = filename(1:LEN_TRIM(filename)),status='old')
    END IF
    
  !use elevation to drift interpolation
  IF ( KeyIsPresent ('elevation-drift', ini, section = 'solar-radiation') ) THEN
	     elevationDrift = IniReadInt ('elevation-drift', ini, section = &
                                      'solar-radiation')
  ELSE
      elevationDrift = 0 !default, suppress drift
  END IF
  
  IF (elevationDrift == 1) THEN
      
     
      IF (interpolationMethod == 0 ) THEN
          CALL Catch ('error', 'SolarRadiation',   &
                'elevation drift cannot be applied when interpolation = 0 ')
      END IF
      
      IF ( KeyIsPresent ('time-zone', ini, section = 'solar-radiation') ) THEN
	     timezone = IniReadInt ('time-zone', ini, section = &
                                      'solar-radiation')
      ELSE
          CALL Catch ('error', 'SolarRadiation',   &
		       'timezone is missing' )
      END IF

      
       !set azimuth angles to compute shadow and convert to radians
       azimuth = (/0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, &
           180.0, 202.5, 225.0, 247.5, 270.0, 292.5, 315.0, 337.5/)
       
       azimuth = azimuth * degToRad
       
       CALL DeriveSlope  (dem, slope)
                    
       CALL DeriveAspect (dem, aspect)
                    
       CALL NewGrid (shadowGrid, dem)
                    
       ALLOCATE ( viewangle (dem % idim,dem % jdim) )
	
  END IF  
  

  !grid exporting settings
  IF (KeyIsPresent ('export', ini, section = 'solar-radiation')  )  THEN
     export = IniReadInt ('export', ini, section = 'solar-radiation')
  ELSE
     export = 0
  END IF 

  IF (export == 1) THEN
     
     !set export path name
     IF (KeyIsPresent ('export-path', ini, section = 'solar-radiation')  )  THEN
         export_path = IniReadString ('export-path', ini, section = 'solar-radiation')
         !check if path ends with / or \
         IF (  export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '\' .AND. &
               export_path (LEN_TRIM (export_path):LEN_TRIM (export_path)) /= '/' ) THEN
             export_path (LEN_TRIM (export_path)+1:LEN_TRIM (export_path)+1) = '\'
         END IF
       
         IF (INDEX (export_path,'.') == 1) THEN !detected relative path, remove '.'
            export_path = export_path (2:LEN_TRIM(export_path))
            !build absolute path
            export_path = TRIM(CurrentDir() ) // TRIM(export_path)
         END IF

         !check OS convention
        IF (GetOS () == WIN32) THEN
          export_path = ReplaceChar (export_path,'/','\')
        ELSE
          export_path = ReplaceChar (export_path,'\','/')
        END IF
  
         
         !check folder exists
         IF ( .NOT. DirExists (TRIM (export_path) ) ) THEN
             CALL Catch ('info', 'SolarRadiation',   &
                  'creating directory:  ', argument = TRIM(export_path))
             CALL DirNew (export_path)
         END IF
     ELSE
         CALL Catch ('error', 'SolarRadiation',   &
                  'missing export-path')
     END IF 
     
     !starting time
     IF (KeyIsPresent ('export-start', ini, section = 'solar-radiation')  )  THEN
         timeString = IniReadString ('export-start', ini, 'solar-radiation')
         export_start = timeString
     ELSE
         CALL Catch ('error', 'SolarRadiation',   &
                  'missing export-start')
     END IF 
     
     !initialize timeNewExport
     timeNewExport = export_start
     
     !ending time
     IF (KeyIsPresent ('export-stop', ini, section = 'solar-radiation')  )  THEN
         timeString = IniReadString ('export-stop', ini, 'solar-radiation')
         export_stop = timeString
     ELSE
         CALL Catch ('error', 'SolarRadiation',   &
                  'missing export-start')
     END IF 
     
     !export dt
     IF (KeyIsPresent ('export-dt', ini, section = 'solar-radiation')  )  THEN
         export_dt = IniReadInt ('export-dt', ini, section = 'solar-radiation')
     ELSE
         CALL Catch ('error', 'SolarRadiation',   &
                  'missing export-dt')
     END IF 

     !coordinate reference system
     IF (KeyIsPresent ('export-epsg', ini, section = 'solar-radiation')  )  THEN
         export_epsg = IniReadInt ('export-epsg', ini, section = 'solar-radiation')
         
         exportGridMapping = DecodeEPSG (export_epsg)
         IF (exportGridMapping == radiation % grid_mapping) THEN
            needConversion = .FALSE.
            !initialize grid for exporting with CRS as meteo variable
            CALL NewGrid (exportedGrid, radiation)
            CALL NewGrid (exportedGridVar, radiation)
            CALL NewGrid (exportedGridNet, netRadiation)
         ELSE
            needConversion = .TRUE.
            !initialize grid for converting with required CRS
            exportedGrid % grid_mapping = DecodeEPSG (export_epsg)
            exportedGridVar % grid_mapping = DecodeEPSG (export_epsg)
            exportedGridNet % grid_mapping = DecodeEPSG (export_epsg)
         END IF
     ELSE
         CALL Catch ('info', 'SolarRadiation',   &
                  'export-epsg not defined, use default')
         needConversion = .FALSE.
         !initialize grid for exporting with CRS 
         CALL NewGrid (exportedGrid, radiation)
         CALL NewGrid (exportedGridVar, radiation)
         CALL NewGrid (exportedGridNet, netRadiation)
     END IF 

     !export file format 
     IF (KeyIsPresent ('export-format', ini, section = 'solar-radiation')  )  THEN
         export_format = IniReadInt ('export-format', ini, section = 'solar-radiation')
     ELSE
         CALL Catch ('error', 'SolarRadiation',   &
                  'missing export-format')
     END IF 

     IF (export_format == 3) THEN
       !grid map  
       CALL SetLongName ( 'downwelling_shortwave_flux_in_air', exportedGrid)
       CALL SetLongName ( 'net_downward_flux_in_air', exportedGridNet)
       CALL SetStandardName ( 'downwelling_shortwave_flux_in_air', exportedGrid)
       CALL SetStandardName ( 'net_downward_flux_in_air', exportedGridNet)
       CALL SetUnits ('Wm-2', exportedGrid) 
       CALL SetUnits ('Wm-2', exportedGridNet) 
       !if file exists, remove it
       export_file = TRIM(export_path) //  'radiation.nc'
       IF ( FileExists (export_file) ) THEN
          CALL FileDelete (export_file)
       END IF
       export_file_net = TRIM(export_path) //  'net_radiation.nc'
       IF ( FileExists (export_file_net) ) THEN
          CALL FileDelete (export_file_net)
       END IF
       
       !variance of kriging 
       CALL SetLongName ( 'downwelling_shortwave_flux_in_air_variance', exportedGridVar)
       CALL SetStandardName ( 'downwelling_shortwave_flux_in_air_variance', exportedGridVar)
       CALL SetUnits ('Wm-2', exportedGridVar) !this unit is for exporting grid
       !if file exists, remove it
       export_file_var = TRIM(export_path) //  'radiation_variance.nc'
       IF ( FileExists (export_file_var) ) THEN
          CALL FileDelete (export_file_var)
       END IF
       
     END IF

  END IF
	
  !complete initialization

  IF (interpolationMethod == 0) THEN
        !Get the dt of imported  field. Assume dt is regular	
        dtGrid = GetDtGrid (filename = fileGrid, checkRegular = .TRUE.)
        !check dt is multiple of dtGrid
        IF (.NOT.(MOD(dtRadiation,dtGrid) == 0)) THEN
            CALL Catch ('error', 'SolarRadiation',   &
                'dt solar radiation is not multiple of dt of input grid')
        END IF
  ELSE
        !populate  metadata
        CALL ReadMetadata (fileunit, radiometers)

        !check spatial reference system
        IF ( .NOT. radiometers % mapping == mask % grid_mapping)  THEN
            CALL Catch ('info', 'SolarRadiation',   &
		              'converting coordinate of stations')
            !convert stations' coordinate
            point1 % system = radiometers % mapping
            point2 % system = mask % grid_mapping
            radiometers % mapping = mask % grid_mapping
            DO i = 1, radiometers % countObs
              point1 % easting = radiometers % obs (i) % xyz % easting
              point1 % northing = radiometers % obs (i) % xyz % northing
              point1 % elevation = radiometers % obs (i) % z
              CALL Convert (point1, point2)
              radiometers % obs (i) % xyz % easting = point2 % easting
              radiometers % obs (i) % xyz % northing = point2 % northing 
            END DO
        END IF
        
  END IF !interpolationMethod
   
RETURN
  END SUBROUTINE SolarRadiationInit

!==============================================================================
!| Description:
!   Read radiation data
!
! References:
!
!  Ranzi R., Rosso R., Un modello idrologico distribuito, su base 
!   fisica, dello scioglimento nivale, Master thesis (in italian), 
!   Politecnico di Milano,1989.
!
SUBROUTINE SolarRadiationRead &
!
( time, dem, albedo, temp, relHum )

IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(IN) :: time !!current time
TYPE (grid_real), INTENT(IN) :: dem !!used to apply drift of station data
TYPE (grid_real), INTENT(IN) :: albedo !!used to apply drift of radiation site data
TYPE (grid_real), INTENT(IN) :: temp !! air temperature (degree celsius)
TYPE (grid_real), INTENT(IN) :: relHum !! air relative humidity (0-100)

!local declarations:
TYPE (DateTime) :: time_toread !! time to start reading from
TYPE (DateTime) :: timeLocal !! local time for applying topographic drift
CHARACTER (LEN = 300) :: filename
CHARACTER (LEN = 300) :: varname
REAL (KIND = float)   :: rsquare
INTEGER (KIND = short) :: i, j
REAL (KIND = float) :: d !!solar declination
REAL (KIND = float) :: w !!solar hour angle
REAL (KIND = float) :: a1 !!sun elevation angle
REAL (KIND = float) :: az !!azimuth
REAL (KIND = float) :: mh !!atmosphere optical depth
REAL (KIND = float) :: z !!terrain elevation
REAL (KIND = float) :: Ic !!clear sky direct radiation
REAL (KIND = float) :: D1 !!clear sky diffuse radiation
REAL (KIND = float) :: DF  !!diffuse radiation from surrounding elements
REAL (KIND = float) :: A  !!reflexed radiation from surrounding elements
REAL (KIND = float) :: Rstar !! clear sky radiation = Ic + D1
REAL (KIND = float) :: radMeas !!measured value of radiation [W/m2]
REAL (KIND = float) :: kt !! clearness index, cloudiness index complement
REAL (KIND = float) :: fB1 !! hillslope coefficient
REAL (KIND = float) :: cosT !! topographic factor
REAL (KIND = float) :: Q !!direct radiation component on inclined surface
REAL (KIND = float) :: radG !! ground radiation

REAL (KIND = float), PARAMETER :: Isc = 1367.  !! solar constant
REAL (KIND = float), PARAMETER :: perclo = 0.22 !!minimum fraction of diffuse radiation
REAL (KIND = float), PARAMETER :: k1 = 0.4  !!atmosphere diffusion coefficient

!-------------------------end of declarations----------------------------------

IF ( .NOT. (time < timeNew) ) THEN
   
   !time_toread = time + - (dtRadiation - radiometers % timeIncrement)
   time_toread = time
   timeString = time_toread
   CALL Catch ('info', 'SolarRadiation',   &
		                 'read new solar radiation data: ', &
                     argument = timeString)

   SELECT CASE (interpolationMethod)
    CASE (0) !input grid
       
	  CALL ReadField (fileGrid,  time_toread, &
         dtRadiation, dtGrid, &
         'M', radiation, &
         varName = radiation % var_name)
          
    CASE DEFAULT !requires interpolation
      !read new  data
      CALL ReadData (network = radiometers, fileunit = fileunit, &
                      time = time_toread, aggr_time = dtRadiation, &
                      aggr_type = 'ave', tresh = valid_prcn)
      
      
      IF (interpolationMethod_assignment == 1) THEN !only one method is applied
                
          CALL Interpolate (network = radiometers, &
                    method = interpolationMethod, &
                    near = neighbors, &
                    idw_power = idw_power, & 
                    anisotropy = krige_anisotropy, &
                    varmodel = krige_varmodel, &
                    lags = krige_lags, &
                    maxlag = krige_maxlag, &
                    grid = radiation, &
                    devst = grid_devst) 
                
      ELSE
          !loop trough interpolation methods
          DO i = 1, 3
             IF (interpolationMethod_vector(i) > 0) THEN
                      
                CALL Interpolate (network = radiometers, &
                        method = interpolationMethod_vector(i), &
                        near = neighbors, &
                        idw_power = idw_power, & 
                        anisotropy = krige_anisotropy, &
                        varmodel = krige_varmodel, &
                        lags = krige_lags, &
                        maxlag = krige_maxlag, &
                        grid =  interpolatedMap (i), &
                        devst = grid_devst) 
                    
             END IF
          END DO

        !compose the final interpolated field
          DO i = 1, interpolationMethod_map % idim
             DO j = 1, interpolationMethod_map % jdim
                IF (interpolationMethod_map % mat(i,j) /= &
                    interpolationMethod_map % nodata ) THEN
                    radiation % mat (i,j) = &
                    interpolatedMap (interpolationMethod_map % mat(i,j)) % mat(i,j)
                END IF
             END DO
          END DO

      END IF  !1 or more interpolation methods   
    
    END SELECT
      
      IF (elevationDrift == 1) THEN !adjust downwelling shortwave 
                                    !radiation to account for topography
          
          timeLocal = time + INT (timezone * hour)
          
          !compute solar declination
          d = SolarDeclination (timeLocal)
          
          !compute solar hour angle
          w = SolarHourAngle (timeLocal)
          
          !compute sun elevation angle
          a1 = SunElevationAngle (timeLocal, latCentroid)
          
          !compute azimuth
          az = SolarAzimuth (timeLocal, latCentroid)
          
          !compute shadow map
          CALL CastShadow (az, a1, viewangle, shadowGrid)
          
          
          DO i = 1, radiation % idim
              DO j = 1, radiation % jdim
                  IF (radiation % mat (i,j) /= radiation % nodata ) THEN
                      
                      radMeas = radiation % mat (i,j)
                      
                      IF ( a1 <= 0 ) THEN !sun is below horizon
                         radiation % mat (i,j) = radMeas
                         netRadiation % mat (i,j) = radMeas
                         netRadiationFAO % mat (i,j) = radMeas
                         CYCLE
                      END IF
          
                      
                      !elevation
                      z = dem % mat (i,j)
                      
                      !optical depth
                      mh = OpticalDepth (timeLocal, latCentroid, z)
                      
                      !clear sky direct radiation
                      Ic = Isc * EXP ( -mh / SIN (a1) ) * SIN (a1)
                      
                      !clear sky diffuse radiation
                      D1 = k1 * ( Isc * SIN (a1) - Ic )
                      
                      !total clear sky radiation
                      Rstar = Ic + D1
                      
                      !clearness index, cloudiness factor complement (Ranzi, 1989)
                      ! kt = 0 fully cloudy, kt = 1 clear sky
                      kt = (  Rstar - radMeas ) / ( ( 1 - perclo ) * Rstar )
		              IF ( kt < 0. ) THEN
			            kt = 0.
		              END IF
		              IF ( kt > 1. ) THEN
			            kt = 1.
                      END IF
                      
                      !diffuse radiation from surrounding elements
		              IF ( Ic < 0. ) THEN
			             DF = 0.
		              ELSE
			             DF = MIN (radMeas, D1 * ( 1 - Kt ) + radMeas * Kt)
                      END IF
                      
                      !reflexed radiation from surrounding elements
                      fB1 = 1. - slope % mat (i,j)
                      A = radMeas * albedo % mat (i,j) * ( 1. - fB1 )
                      
                      !direct radiation component
                      cosT = COS (a1) * SIN (slope % mat(i,j) ) * &
                             COS ( az - aspect % mat(i,j) ) + &
                             SIN (a1) * COS (slope % mat(i,j) )
                      
                      Q = ( radMeas - DF ) * cosT / SIN (a1)
                      
                      !check shadow
                      IF ( shadowGrid % mat(i,j) == 1 .OR. Q < 0 ) THEN
			                radG = DF + A
			                cosT = -1
                      !error computing radG
                      ELSE IF( ( Q + DF + A ) > Isc) THEN
			                radG = Isc
                      ELSE
			                radG = Q + DF + A
                      END IF
                        
                      radiation % mat (i,j) = radG     
                      
                      !compute net radiation
                      CALL ComputeNetRadiation ( kt, albedo % mat (i,j), &
                                 radiation % mat (i,j), temp % mat(i,j),  &
                                 relHum % mat (i,j), netRadiation % mat (i,j) )
                      
                  END IF
              END DO
          END DO
          
      ELSE !compute only cloudiness factor and net radiation
          
          timeLocal = time + INT (timezone * hour)
          
          !compute sun elevation angle
          a1 = SunElevationAngle (timeLocal, latCentroid)
          
          DO i = 1, radiation % idim
              DO j = 1, radiation % jdim
                  IF (radiation % mat (i,j) /= radiation % nodata ) THEN
                      
                      radMeas = radiation % mat (i,j)
                      
                      IF ( a1 <= 0 ) THEN !sun is below horizon
                         netRadiation % mat (i,j) = radMeas
                         netRadiationFAO % mat (i,j) = radMeas
                         CYCLE
                      END IF
                      
                      !elevation
                      z = dem % mat (i,j)
                      
                      !optical depth
                      mh = OpticalDepth (timeLocal, latCentroid, z)
                      
                      !clear sky direct radiation
                      Ic = Isc * EXP ( -mh / SIN (a1) ) * SIN (a1)
                      
                      !clear sky diffuse radiation
                      D1 = k1 * ( Isc * SIN (a1) - Ic )
                      
                      !total clear sky radiation
                      Rstar = Ic + D1
                      
                      !cloudiness factor complement.
                      ! kt = 0 fully cloudy, kt = 1 clear sky
                      kt = (  Rstar - radMeas ) / ( ( 1 - perclo ) * Rstar )
		              IF ( kt < 0. ) THEN
			            kt = 0.
		              END IF
		              IF ( kt > 1. ) THEN
			            kt = 1.
                      END IF
                      
                      !compute net radiation
                      CALL ComputeNetRadiation ( kt, albedo % mat (i,j), &
                                 radMeas, temp % mat(i,j), relHum % mat (i,j), &
                                 netRadiation % mat (i,j) )
                      
                      !compute net radiation FAO with albedo = 0.23
                      CALL ComputeNetRadiation ( kt, 0.23, &
                                 radMeas, temp % mat(i,j), relHum % mat (i,j), &
                                 netRadiationFAO % mat (i,j) )
        
                      
                  END IF
              END DO
          END DO
          
      END IF  !elevationDrift
		
	! END SELECT moved up

  !apply scale factor and offset, if defined
	IF (offset_value /= MISSING_DEF_REAL) THEN
	   CALL Offset (radiation, offset_value)
       CALL Offset (netRadiation, offset_value)
	END IF
	
	
	IF (scale_factor /= MISSING_DEF_REAL) THEN
	   CALL Scale (radiation, scale_factor)
       CALL Scale (netRadiation, scale_factor)
    END IF


  !grid exporting
  IF(export > 0 ) THEN
	  IF( time == timeNewExport .AND. &
        time >= export_start .AND. &
        time <= export_stop ) THEN
        timeString = time
        timeString = timeString(1:19)
        timeString(14:14) = '-'
		    timeString(17:17) = '-'
        
            
        !radiation
        grid_devst % reference_time = radiation % reference_time
        IF (needConversion) THEN
           CALL GridConvert (radiation, exportedGrid)
           CALL GridConvert (grid_devst, exportedGridVar)
        ELSE
           exportedGrid = radiation
           exportedGridVar = grid_devst
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_radiation.asc'
              CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_ASCII)
              
              IF (krige_var == 1) THEN !export kriging variance
                    export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                           '_radiation_variance.asc'
                    CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                    CALL ExportGrid (exportedGridVar, export_file_var, ESRI_ASCII)
              END IF
              
        CASE (2) !esri-binary
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_radiation'
              CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, ESRI_BINARY)
              
              IF (krige_var == 1) THEN !export kriging variance
                   export_file_var = TRIM(export_path) //  TRIM(timeString) // &
                               '_radiation_variance'
                   CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting variance grid to file: ', &
                              argument = TRIM(export_file_var))
                   CALL ExportGrid (exportedGridVar, export_file_var, ESRI_BINARY)
              END IF
              
        CASE (3) !net_cdf
              CALL SetCurrentTime (time, exportedGrid)
              CALL Catch ('info', 'SolarRadiation',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGrid, export_file, 'radiation', 'append')
              
              IF (krige_var == 1) THEN !export kriging variance
                  CALL SetCurrentTime (time, exportedGridVar)
                  CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting grid to file: ', &
                              argument = TRIM(export_file_var))
                 CALL ExportGrid (exportedGridVar, export_file_var, &
                                  'radiation_variance', 'append')
              END IF
        END SELECT
        
        
        
        !net radiation
        IF (needConversion) THEN
           CALL GridConvert (netRadiation, exportedGridNet)
        ELSE
           exportedGridNet = netRadiation
        END IF 

        SELECT CASE (export_format)
        CASE (1) !esri-ascii
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_net_radiation.asc'
              CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGridNet, export_file, ESRI_ASCII)
              
              
        CASE (2) !esri-binary
              export_file = TRIM(export_path) //  TRIM(timeString) // &
                           '_net_radiation'
              CALL Catch ('info', 'SolarRadiation',   &
		                       'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGridNet, export_file, ESRI_BINARY)
              
              
        CASE (3) !net_cdf
              CALL SetCurrentTime (time, exportedGridNet)
              CALL Catch ('info', 'SolarRadiation',   &
                           'exporting grid to file: ', &
                           argument = TRIM(export_file))
              CALL ExportGrid (exportedGridNet, export_file_net, 'net_radiation', 'append')
              
              
        END SELECT
        
        
        timeNewExport = timeNewExport + export_dt
    END IF
  ENDIF

  !time forward
  timeNew = timeNew + dtRadiation
END IF

RETURN
END SUBROUTINE SolarRadiationRead


!==============================================================================
!| Description:
!   set properties and initialize variables for each interpolation method
SUBROUTINE SetSpecificProperties &
!
( method, ini, mask )

IMPLICIT NONE

!Arguments with intent(in):
INTEGER (KIND = short), INTENT(in) :: method
TYPE (IniList)        , INTENT(in) :: ini
TYPE (grid_integer)   , INTENT(in) :: mask

!----------------------------end of declarations-------------------------------

!set data specific for each interpolation method
	SELECT CASE (method)
    CASE (1) !thiessen
      CALL NewGrid (interpolatedMap (1), mask, 0.)
    CASE (2) !inverse distance
      CALL NewGrid (interpolatedMap (2), mask, 0.)
	    !nearest-points
      IF (KeyIsPresent ('nearest-points', ini, &
                         section = 'solar-radiation') ) THEN
	       neighbors = IniReadInt ('nearest-points', ini, &
                                   section = 'solar-radiation')
      ELSE !nearest-points missing
          CALL Catch ('error', 'SolarRadiation',   &
		                 'nearest-points missing in meteo configuration file')
      END IF
 
    CASE (3) !kriging
      CALL NewGrid (interpolatedMap (3), mask, 0.)
      !nearest-points
      IF (KeyIsPresent ('nearest-points', ini, &
                         section = 'solar-radiation') ) THEN
	       neighbors = IniReadInt ('nearest-points', ini, &
                                   section = 'solar-radiation')
           CALL Catch ('info', 'SolarRadiation', 'neighbors set to: ', &
                        argument = ToString (neighbors) )
      ELSE !nearest-points missing
          CALL Catch ('error', 'SolarRadiation',   &
		                 'nearest-points missing in meteo configuration file')
      END IF
      
      !kriging-variance
      IF (KeyIsPresent ('kriging-variance', ini, &
                         section = 'solar-radiation') ) THEN
	       krige_var = IniReadInt ('kriging-variance', ini, &
                                   section = 'solar-radiation')
           
           CALL Catch ('info', 'SolarRadiation',   &
		                 'kriging variance export set to ', &
                          argument = ToString(krige_var) )
      ELSE !kriging-variance missing
           krige_var = 0                                     
           CALL Catch ('info', 'SolarRadiation',   &
		                 'kriging variance export set to default (0)')
      END IF
      
      !kriging-anisotropy
      IF (KeyIsPresent ('kriging-anisotropy', ini, &
                         section = 'solar-radiation') ) THEN
	       krige_anisotropy = IniReadInt ('kriging-anisotropy', ini, &
                                   section = 'solar-radiation')
           CALL Catch ('info', 'SolarRadiation', 'kriging anisotropy set to: ', &
                       argument = ToString (krige_anisotropy) )
      ELSE !kriging-anisotropy missing
           krige_anisotropy = 0                                     
           CALL Catch ('info', 'SolarRadiation',   &
		                 'kriging anisotropy set to default (0)')
      END IF
      
      !kriging-variogram-model
      IF (KeyIsPresent ('kriging-variogram-model', ini, &
                         section = 'solar-radiation') ) THEN
	       krige_varmodel = IniReadInt ('kriging-variogram-model', ini, &
                                   section = 'solar-radiation')
           CALL Catch ('info', 'SolarRadiation', 'kriging variogram model set to: ', &
                       argument = ToString (krige_varmodel) )
           IF (krige_varmodel == 0 ) krige_varmodel = 5 !automatic fitting
      ELSE !kriging-variogram-model
           krige_varmodel = 2                                     
           CALL Catch ('info', 'SolarRadiation',   &
		                 'kriging variogram model set to default (2 exponential)')
      END IF
      
      !kriging-lags
      IF (KeyIsPresent ('kriging-lags', ini, &
                         section = 'solar-radiation') ) THEN
	       krige_lags = IniReadInt ('kriging-lags', ini, &
                                   section = 'solar-radiation')
           IF (krige_lags == 0) krige_lags = 15
           CALL Catch ('info', 'SolarRadiation', 'kriging lags set to: ', &
                       argument = ToString (krige_lags) )
      ELSE !kriging-variogram-model
           krige_lags = 15                                     
           CALL Catch ('info', 'SolarRadiation',   &
		                 'kriging lags set to default (15)')
      END IF
      
       !kriging-maxlag
      IF (KeyIsPresent ('kriging-maxlag', ini, &
                         section = 'solar-radiation') ) THEN
	       krige_maxlag = IniReadInt ('kriging-maxlag', ini, &
                                   section = 'solar-radiation')
           CALL Catch ('info', 'SolarRadiation', 'kriging max lag set to: ', &
                     argument = ToString (krige_maxlag) )
      ELSE !kriging-maxlag missing
           krige_maxlag = 0
          CALL Catch ('error', 'SolarRadiation',   &
		                 'kriging max lag set to default (0)')
      END IF
	
	END SELECT


RETURN
END SUBROUTINE SetSpecificProperties


!==============================================================================
!| Description:
!   Compute the maximum angle of sky obstruction along 16 directions
!
! References:
!
! Zhang, Y., Chang, X. & Liang, J. Comparison of different algorithms for 
!  calculating the shading effects of topography on solar irradiance in 
!  a mountainous area. Environ Earth Sci 76, 295 (2017). 
!  https://doi.org/10.1007/s12665-017-6618-5
!   
SUBROUTINE SkyView &
!
( dem, azimuth, view)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (grid_real)   , INTENT(in) :: dem
REAL (KIND = float), INTENT(in) :: azimuth (:)  !azimuth angles [rad]

!Arguments with intent(inout):
TYPE (ViewingAngle), INTENT(inout) :: view(:,:)

!local declarations
INTEGER (KIND = short) :: i,j,l,i1,j1
TYPE (Coordinate) :: x0y0, x1y1
REAL (KIND = float) :: delta, deltax, deltay
REAL (KIND = float) :: length
LOGICAL :: check
REAL (KIND = float) :: topAngle, topAngleMax

!----------------------------end of declarations-------------------------------
x0y0 % system = dem % grid_mapping
x1y1 % system = dem % grid_mapping
delta = dem % cellsize / 2.

DO i = 1, dem % idim
    DO j = 1, dem % jdim
        IF (dem % mat (i,j) /= dem % nodata) THEN
            DO l = 1, SIZE (azimuth)
                !get coordinate of starting point
                CALL GetXY(i, j, dem, x0y0 % easting, x0y0 % northing)
        
                deltax = 0.
                deltay = 0.
                topAngle = 0.
                topAngleMax = 0.
                DO 
                  !add increment in x and y
                  deltax = deltax + delta * SIN (azimuth (l) )
                  deltay = deltay + delta * COS (azimuth (l) )
                  x1y1 % easting = x0y0 % easting + deltax
                  x1y1 % northing = x0y0 % northing + deltay
                    
                  !retrieve local coordinate of new cell
                  CALL GetIJ (x1y1 % easting, x1y1 % northing, &
                              dem, i1, j1, check)
                  
                  IF (.NOT. check ) THEN
                      EXIT
                  ELSE IF (dem % mat (i1, j1) == dem % nodata ) THEN
                      EXIT
                  END IF
                  
                  !compute distance 
                  length = Distance (x1y1, x0y0)
                  
                  !compute topographic angle
                  topAngle = ATAN (  (dem % mat (i1, j1) - dem % mat (i,j) ) &
                                     / length )
                  IF (topAngle > topAngleMax) topAngleMax = topAngle
                END DO
                view(i,j) % angle (l) = topAngleMax
            END DO
        END IF
    END DO
END DO

RETURN
END SUBROUTINE SkyView



!==============================================================================
!| Description:
!   Compute shadow grid
SUBROUTINE CastShadow &
!
( az, sunHeight, view, grid)
    
IMPLICIT NONE

!Arguments with intent(in):

REAL (KIND = float), INTENT(in) :: az !current sun azimuth [rad]
REAL (KIND = float), INTENT(in) :: sunHeight !current sun height [rad]
TYPE (ViewingAngle), INTENT(in) :: view (:,:)

!Arguments with intent(inout):
TYPE (grid_real)   , INTENT(inout) :: grid


!local declarations
INTEGER (KIND = short) :: i, j, m, l
REAL (KIND = float) :: azDiff, azDiffMin

!----------------------------end of declarations-------------------------------

!search for closer azimuth to current sun azimuth
azDiffMin = 100.
l = 0
DO m = 1, SIZE (azimuth)
    azDiff = ABS ( azimuth (m) - az)
    IF ( azDiff < azDiffMin) THEN
        azDiffMin = azDiff
        l = m
    END IF
END DO

DO i = 1, grid % idim
    DO j = 1, grid % jdim
        IF ( grid % mat (i,j) /= grid % nodata ) THEN
            IF (sunHeight < view (i,j)  % angle (l) ) THEN
                grid % mat (i,j) = 1
            ELSE
                grid % mat (i,j) = 0
            END IF
        END IF
    END DO
END DO

RETURN
END SUBROUTINE CastShadow


!==============================================================================
!| Description:
!   Compute solar declination. The declination of the sun is the angle 
!   between the equator and a line drawn from the centre of the Earth to 
!   the centre of the sun.
! 
! References:
!
!   Iqbal, M.: An introduction to solar radiation, Academis Press
!     Canada, Ontario, 1983.
FUNCTION SolarDeclination &
!
(time) &
!
RESULT (sdec)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(in) :: time

!local declarations:
REAL (KIND = float) :: sdec ![radians]
REAL (KIND = float) :: dayAngle ![radians]
INTEGER (KIND = short) :: dn !day number of the year


!------------------------------------end of declarations-----------------------

!day of year
dn = DayOfYear (time, 'noleap')

!day angle
dayAngle = 2. * pi * (dn -1)  / 365.

!declination
sdec = (0.006918 - 0.399912 * COS (dayAngle) + 0.070257 * SIN (dayAngle) - &
        0.006758 * COS (2 * dayAngle) + 0.000907 * SIN (2 * dayAngle) - &
        0.002697 * COS (3 * dayAngle) + 0.00148 * SIN (3 * dayAngle))
RETURN
END FUNCTION SolarDeclination


!==============================================================================
!| Description:
!   Compute the sun elevation angle
!
! References:
!
!   Gates DM (1980) Biophysical ecology. Springer, New York,
!   page 101, eq. 6.6
FUNCTION SunElevationAngle &
!
(time, lat) &
!
RESULT (sea)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(in) :: time
REAL (KIND = float), INTENT(in) :: lat !latitude [radians]

!local declarations:
REAL (KIND = float) :: sea ! sun elevation angle [radians]
REAL (KIND = float) :: dec !solar declination [radians]
REAL (KIND = float) :: w !hour angle [radians]
!------------------------------------end of declarations-----------------------

w = SolarHourAngle (time)

dec = SolarDeclination (time)
sea = ASIN ( COS (lat) * COS (dec) * COS (w) + SIN (lat) * SIN (dec) )

RETURN
END FUNCTION SunElevationAngle


!==============================================================================
!| Description:
!  Compute the hour angle.  the solar hour angle is an expression of time, 
!   hour angle is 0.000 degree, with the time before solar noon
!   expressed as negative degrees, and the local time after solar noon
!   expressed as positive degrees. For example, at 10:30 AM local
!   apparent time the hour angle is -22.5 degree
!   (15 degree per hour times 1.5 hours before noon).
!
! References:
!
!   Iqbal, M.: An introduction to solar radiation, Academis Press
!     Canada, Ontario, 1983.
FUNCTION SolarHourAngle &
!
(time) &
!
RESULT (w)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(in) :: time

!local declarations:
REAL (KIND = float) :: w ! hour angle [radians]
INTEGER (KIND = short) :: t !hour
!------------------------------------end of declarations-----------------------

t = GetHour (time)
w = 15. * pi / 180. * (t - 12.)

RETURN
END FUNCTION SolarHourAngle


!==============================================================================
!| Description:
!  Compute tthe atmospheric optical depth
!
! References:
!
!   Kreith, F., and Kreider, J. F., 1978, Principles of Solar Engineering 
!     (New York: McGraw-Hill ).
!
!   Cartwright, T. J., 1993, Modelling the World in a Spreadsheet: 
!      Environmental Simulation on a Microcomputer 
!      (Baltimore: Johns Hopkins University Press).
!
!   LALIT KUMAR, ANDREW K. SKIDMORE & EDMUND KNOWLES (1997) Modelling 
!     topographic variation in solar radiation in a GIS environment, 
!     International Journal of Geographical Information Science, 
!     11 :5, 475-497, DOI: 10.1080/136588197242266
FUNCTION OpticalDepth &
!
(time, lat, z) &
!
RESULT (s)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(in) :: time
REAL (KIND = float), INTENT(in) :: lat !latitude [radians]
REAL (KIND = float), INTENT(in) :: z !terrain elevation [m]

!local declarations:
REAL (KIND = float) :: s !!optical depth 
REAL (KIND = float) :: presscorr !atmospheric correction factor
REAL (KIND = float) :: s0 !depth at sea level
REAL (KIND = float) :: a1 ! sun elevation angle [radians]
!------------------------------------end of declarations-----------------------

!compute sun elevation angle
a1 = SunElevationAngle (time, lat)


!compute atmospheric optical depth at sea level
s0 = ( 1229. + (614. * SIN (a1) ) **2. )**0.5 -614. * SIN (a1)

!compute pressure correction factor
presscorr = ( ( 288.15 - 0.0065 * z ) / 288.15 ) ** 5.25588

!optical depth
s = s0 * presscorr

RETURN
END FUNCTION OpticalDepth


!==============================================================================
!| Description:
!  Compute azimuth angle of the Sun's position in the north-clockwise 
!  convention
!
! References:
!
!   Oke, T.R., Boundary layer climates, Second edition, Routledge, 1987.
!    Appendix A1, eq. A1.2
FUNCTION SolarAzimuth &
!
(time, lat) &
!
RESULT (az)
    
IMPLICIT NONE

!Arguments with intent(in):
TYPE (DateTime), INTENT(in) :: time
REAL (KIND = float), INTENT(in) :: lat !latitude [radians]

!local declarations:
REAL (KIND = float) :: az !solar azimuth [radians]
REAL (KIND = float) :: w ! hour angle [radians]
REAL (KIND = float) :: d ! solar declination [radians]
REAL (KIND = float) :: a1 ! sun elevation angle [radians]
INTEGER (KIND = short) :: t !hour
REAL (KIND = float) :: term 
!------------------------------------end of declarations-----------------------

!compute hour angle
w = SolarHourAngle (time)

!compute declination
d = SolarDeclination (time)

!compute sun elevation angle
a1 = SunElevationAngle (time, lat)


!get hour
t = GetHour (time)

term = ( SIN (d) * COS (lat) - COS (d) * SIN (lat) * COS (w) ) / COS (a1)

! compute azimuth

IF ( term < -1.) THEN
    az = pi
ELSE IF (term > 1.) THEN
    az = 0.
ELSE IF (t <= 12) THEN
    az = ACOS ( term )
ELSE
    az = 2. * pi - ACOS ( term )
END IF
    
RETURN
END FUNCTION SolarAzimuth


!==============================================================================
!| Description:
!  Compute net radiation
!
! References:
!
!  Wales-Smith, B. G. (1980). Estimates of net radiation for evaporation 
!   calculations, Hydrological Sciences Journal, 25:3, 237-242, 
!   DOI: 10.1080/02626668009491931
!
!   
SUBROUTINE ComputeNetRadiation &
!
(cfc, alb, short, hum, temp, netRad)

    
IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(in) :: cfc !!cloudiness factor complement 
REAL (KIND = float), INTENT(in) :: alb !!albedo
REAL (KIND = float), INTENT(in) :: short !!shortwave radiation (W/m2)
REAL (KIND = float), INTENT(in) :: hum !! air relative hunidity (0-100)
REAL (KIND = float), INTENT(in) :: temp !! air temperature (degree celsius)

!Arguments with intent(out):
REAL (KIND = float), INTENT(out) :: netRad !!net radiation


!local declarations:
REAL (KIND = float) :: cf !!cloudiness factor
REAL (KIND = float) :: netShort !!net shortwave radiation (W/m2)
REAL (KIND = float) :: netLong !!net longwave radiation (W/m2)
REAL (KIND = float) :: albedotemp !!temporary albedo value
REAL (KIND = float) :: hum01 !!air relative humidity (0-1)
REAL (KIND = float) :: humMin !!minimum value for air relative humidity (0-100)
REAL (KIND = float) :: es !! saturation vapor pressure (Pa)
REAL (KIND = float) :: ea !!actual vapor pressure (Pa)
!------------------------------------end of declarations-----------------------

!cloudiness factor
cf = 1. - cfc

!net shortwave radiation
IF ( alb < 0. ) THEN
	albedotemp = 0.2
	netShort = (1. - albedotemp) * short
ELSE
	netShort = ( 1. - alb ) * short
END IF

!vapor pressure
humMin = 10
IF ( hum > humMin ) THEN
    hum01 = hum / 100.
ELSE
    hum01 = humMin / 100.
END IF

es = 6.107 * exp ( ( 17.27 * temp ) / ( temp + 237.3 ) )
ea = hum01 * es

!net longwave radiation ( Wales-Smith, 1980)

netLong = - ( stefanBoltzman * (temp + 273.15 )**4. ) * &
           ( 0.56 - 0.079 * ea**0.5 ) * ( 0.1 + 0.9 * cf )
	
!total net radiation
netRad = netShort + netLong

RETURN
END SUBROUTINE ComputeNetRadiation


END MODULE SolarRadiation