SolarRadiationRead

public subroutine SolarRadiationRead(time, dem, albedo, temp, relHum)

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.

Arguments

Type	Intent		Name
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)

Variables

Type	Visibility	Attributes		Name		Initial
real(kind=float),	public		::	A			reflexed radiation from surrounding elements
real(kind=float),	public		::	D1			clear sky diffuse radiation
real(kind=float),	public		::	DF			diffuse radiation from surrounding elements
real(kind=float),	public		::	Ic			clear sky direct radiation
real(kind=float),	public,	parameter	::	Isc	=	1367.	solar constant
real(kind=float),	public		::	Q			direct radiation component on inclined surface
real(kind=float),	public		::	Rstar			clear sky radiation = Ic + D1
real(kind=float),	public		::	a1			sun elevation angle
real(kind=float),	public		::	az			azimuth
real(kind=float),	public		::	cosT			topographic factor
real(kind=float),	public		::	d			solar declination
real(kind=float),	public		::	fB1			hillslope coefficient
character(len=300),	public		::	filename
integer(kind=short),	public		::	i
integer(kind=short),	public		::	j
real(kind=float),	public,	parameter	::	k1	=	0.4	atmosphere diffusion coefficient
real(kind=float),	public		::	kt			clearness index, cloudiness index complement
real(kind=float),	public		::	mh			atmosphere optical depth
real(kind=float),	public,	parameter	::	perclo	=	0.22	minimum fraction of diffuse radiation
real(kind=float),	public		::	radG			ground radiation
real(kind=float),	public		::	radMeas			measured value of radiation [W/m2]
real(kind=float),	public		::	rsquare
type(DateTime),	public		::	timeLocal			local time for applying topographic drift
type(DateTime),	public		::	time_toread			time to start reading from
character(len=300),	public		::	varname
real(kind=float),	public		::	w			solar hour angle
real(kind=float),	public		::	z			terrain elevation

Source Code

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

SolarRadiationRead Subroutine

Contents

Variables

Source Code

public subroutine SolarRadiationRead(time, dem, albedo, temp, relHum)

Arguments

Variables

Source Code