[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[TRNSYS-users] Type109 high latitude bug



Hi Fan and TRNSYS users,

Sorry for the late reply, it's been midsummer holidays...

A bug was discovered in Type109 just three weeks ago. Locations above
approximately latitude 63 deg will have a radiation blackout during
summer months. Apparently there was some sunset checking that caused an
error with locations with midnight sun during summer.

Diego fixed the bug, but the fix is not included in release 16.01.0003
as it was discovered after the release.
I include a fixed Type109 source code, which Diego sent me. You have to
rebuild TRNdll.dll with this replaced Type109 source code file.

Another solution is that you use Type15-2 which does not have this bug.


Kind regards,
Janne

-------------------------------------------------------------
Janne Paavilainen MSc, Doktorand
Forskare inom Energi och miljöteknik
-------------------------------------------------------------
Solar Energy Research Center SERC
Högskolan Dalarna
SE-781 88 Borlänge
-------------------------------------------------------------
tel +46 (0)23 778728
fax +46 (0)23 778701
e-mail: jip@du.se
-------------------------------------------------------------
www.du.se
www.serc.se
www.eses.org
www.uni-kassel.de/fb15/ite/solar/solnet/
-------------------------------------------------------------




Jianhua Fan wrote:

Dear Diego,


I have posted this help for some days but without answer. Can you have a
look at my problem and point the way out for me?

The problem is with type 109 user format. The 109 user type works all
right with my weather data for Copenhagen with the latitude approx. 56°,
while it is very strange with my weather data for Greenland with the
latitude 70.7°. The problem is that there is no solar radiation in a
period in summer when it is supposed to have solar radiation. The
problem will disappear by simply changing the latitude to 56, but that
is not the right solution.

The Trnsys I am using is 16.01.0002.

I am frustrated with type 109 user and thinking about switching back to
TRN 14 which I am sure does not have this kind of stupid problem.

Herewith I send you all the tpf file and the weather data file in view
that some of you might know the reason. In order to run the simulation,
two nonstandard dlls are needed. They are also enclosed in the email.

I am looking forward to hearing from you.

Thank you in advance.



Best regards

Jianhua Fan

______________________________________

Department of Civil Engineering

Technical University of Denmark

Building 118, Brovej

DK-2800 Kgs.Lyngby

Denmark



Phone: +45 45 25 18 89

Fax: +45 45 93 17 55

E-mail: jif@byg.dtu.dk <mailto:jif@byg.dtu.dk>

______________________________________

! TRNSYS Type 109: Combined Data reader and radiation processor
! ----------------------------------------------------------------------------
!
! This routine reads a data file in user format or one of the standard formats
! recognized by Trnsys and performs solar radiation processing 
!
! SubPrograms
! ---------------------------------------------------------------------------
! 
! T109Data:         Data module
! ProcessRadiation: Process the solar radiation
!*****************************************************************************
! Level1	Level2	Level 3
!*****************************************************************************
! TYPE109
!----------INITIAL TIMESTEP
!      		inip
!			input_init
!			infiletest
!			dataread
!
!----------Each Call
!			update_time_sun_radiation
!			akimaip
!			rad_proc
!					beamdifftotal
!					calc_every_step
!					no_radiation
!					sloped_surfaces
!					sloped_surfaces
!					suncalc
!
! Required external libraries
! ----------------------------------------------------------------------------
! None
!
! ----------------------------------------------------------------------------------------------------------------------
! Copyright © 2006 Solar Energy Laboratory, University of Wisconsin-Madison. All rights reserved.
!
! Revision history
! ---------------------------------------------------------------------------------------------------------------------
!
! This type was originally written by Thomas Lechner at TRANSSOLAR Energietechnik GmbH, from March 1996 to February 2001
!
! Modifications:
! 2007-05-30 - DAA - Removed a check for sunset that was conflicting with the solar radiation in high latitude locations.
! 2006-05-09 - MH  - Missing initialization caused error if no wind data was read: u(m)%ro=0.; u(m)%ipol=0.; u(m)%prodf=0.; u(m)%addf=0.; u(m)%samp=0.) 
!                  - warning  text corrected:  1550  FORMAT
! 2007-01-21 - DAA - Fixed a division by zero under some conditions with zero radiation.
! 2007-01-19 - DAA - Added a case for user-defined weather data in which the radiation inputs are beam normal radiation
!                    and diffuse radiation on horizontal.
! 2006-11-29 - DAA - Allowed for reading the slope and azimuth at each time step, using tracking mode 1.
! 2006-09-11 - DAA - Added a call to update_time_sun_radiation during initialization call, in order to initialize the
!                     sunrise and sunset angles appropriately.
! 2006-07-27 - DAA - Added a check if the axis of the surface is horizontal in TRACKING mode 3.
!                    Undid change on 2005-10-25 that swaped equations for TRACKING modes 2 and 3.
! 2006-06-29 - DAA - u(m)%samp was not initialized correctly in user-specified weather data.  The initialization fixed
!                     a bug in Release mode at the first time step.
! 2006-03-01 - DAA - Fixed a bug at sunset that caused an array to be called with an index equal to zero.
! 2006-01-24 - DAA - Changed the way MESSAGES are handled, from STOP to standard CALL MESSAGES subroutine.  
!                    Deleted ERRORFLAG, since it was never used
! 2005-12-05 - MH  - Bug fix 1129 (different radiation values in the 2nd year due to lack of accuracy)
! 2005-11-29 - MH/THL - bug fix 1109 (Equations of tracking mode 2 and 3 were mixed up)
!                     - day=MOD(ceiling(tims/24.),365) insert MOD because decl differed for the 2nd year due to a lack of accuracy
! 2005-10-25 - YY/DAA - Added the recalling of the previous values of Tambient and RH.
!                    The values of the last interpolation are stored in an array, in order to keep the values for 
!                    different instances of TYPE109.
! 2004-10-28 - MKu - Handle Southern hemisphere latitude values properly in TMY2 mode(were not present in the original 
!                    US-TMY2 stations but happen in Meteonorm-generated files)
!  TL, 02.03.05, singular values in radiation treatment fixed, enhanced speed due to pointers in data-array u%p
!  TL, 10.05.04, two radiation data files did not work, sunsetflag had to be changed to u(m)%sunsetflag
!		 for different calling units u(m)
!  TL, 26.02.03, new output added: pressure in pa, only mode 2/3 for use in ventilation (comis)
!  TL, 26.02.03, 2 x Mode0 with comment line in data file resulted in error, some variables had
!			   to be reset for each distinct unit, fixed
!  TL, 24.02.03, Mode0 starting after beginning of year while rewind option active
!			   results in zero values-fixed
!  TL, 23.02.03, Parameters had not been checked for wrong inputs, fixed,
!			   while debugging, MYSTOP did not work properly: therefore completely sustituted
!			   by FORTRAN-STOP!!!
!  TL, 22.10.02, TRY-Format was not the original - test ok
!  TL, 21.05.02, extreme radiation peak appears with compiler change not solved jet???
!  TL, 30.04.02, some trouble with radiation data at night - test ok
!                renormalisation of radiation when radiation=0 during day time adjusted
!  TL, 15.04.02, renormalisation of radiation data after interpolation adjusted
!  TL, 09.04.02, peaks in radiation in time step of sunrise --> solar radiation based on sim. time step
!  TW 18.03.2002: comment for "icol" replaced row by col
!  TL, 11.03.02, logical errors in allocation of radiation result array fixed
!  TL, 11.03.02, Bugs due to number and sequence of data columns fixed
!  TL, 15.10.01, Bugs related to modes "Reindl" and "reduced Reindl"  fixed
!  TL, 02.10.01, Mode 1 with IGlob_H, IBeam_N as Input --> Problem Id=0 fixed
!  TL, 29.09.01, Sky Radiation mode99 implemented --->  radiation data set to zero
!  TL 03.09.2001: "irow" replaced by "icol"
!  TL, 03.09.01, checks for negative radiation data implemented
!  TL, 02.09.01, revision of interpolation method (when data rewinds, the last and
!                the first data weren't interpolated until now)
!  TL, xx.xx.01, solar time implemented as a mode 1 definition:  <gmt>     solar
!  TL, 02.05.01, rewind not properly working problem fix
!
! ---------------------------------------------------------------------------------------------------------------------
       
module Type109RadiationData

  real        :: cdelt        !1/100 TH OF DELTA USED TO DETERMINE WHEN SUN IS DOWN
  real        :: hdelt        !HALF OF DELTA used in subroutine sol_rad_proc

  type rad_proc_call
    integer   :: designday                !fixed design days for vdi 2078
    integer   :: gmt                      !time shift in hours from GMT, east: positive
    integer   :: upper                    !UPPER AREA BOUND FOR RESULTS
    integer   :: ias                      !treat as solar: time of weather data means solar time
    integer   :: infloc(2)                !unitnumber, type number
    integer   :: ivdi                     !Rechenverfahren nach VDI
    integer   :: nerr                     !number of ERRORS "RADIATION EXCEED"
    integer   :: ninp                     !number of inputs
    integer   :: nsurf                    !number of surfaces
    integer   :: radmode                  !horiz radiation mode for rad proc
    integer   :: skymode                  !sky model for radiation processor
    integer   :: trackmode                !tracking mode for rad proc
    real      :: iges                     !RADIATION ON HORIZONTAL
    real      :: ib                       !BEAM RADIATION ON HORIZONTAL
    real      :: id                       !DIFFUSE RADIATION ON HORIZONTAL
    real      :: ibdn                     !BEAM RADIATION ON DIRECT NORMAL
    real      :: latit                    !latitude used for radiation data calculation, radians
    real      :: rho                      !GROUND REFLECTANCE
    real      :: rh_amb                   !AMBIENT HUMIDITY FOR FULL REINDL, MODE 6
    real      :: shift                    !shift between local time and solar time, degrees
    real      :: sinlat, coslat, tanlat   !ANGLE DEPENDENT QUANTITIES
    real      :: t_amb                    !AMBIENT TEMPERATURE FOR FULL REINDL
    real, dimension (:), pointer ::slopes !SLOPES OF SURFACES
    real, dimension (:), pointer :: axes  !AZIMUTHS OF SURFACES
  end type rad_proc_call

  type readvar                            !HELPS TO READ USER DEFINED DATA READ
    character (len=100)   :: name, pname  !VARIABLE NAME
    integer               :: icol         !WHICH COL CONTAINS VAR?
    integer               :: ipol, samp   !INTERPOLATION MODE, SAMPLING RATE
    real                  :: mult, add    !MODIFICATION FACTORS
  end type readvar

 end module Type109RadiationData
! -----------------------------------------------------------------------------------
! -----------------------------------------------------------------------------------


SUBROUTINE TYPE109(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
! Export this subroutine for its use in external DLLs
!dec$attributes dllexport :: type109

!************************************************************************************
!   THIS ROUTINE READS IN DATA FROM A FILE, TRANSFORMS
!   THEM BY MULTIPLICATION AND ADDITION AND PROVIDES 3 MODES OF INTERPOLATION
!   IF WEATHER DATA INCLUDING RADIATION ARE READ IN, A RADIATION PROCESSOR SUBROUTINE
!   CALCULATES DIRECT AND DIFFUSE RADIATION INCIDENCE ON AN ARBITRARY NUMBER OF
!   SURFACES.
!
!   NOTE THAT ONLY MAXU TYPE 109 UNITS ARE ALLOWED
!
!   THIS ROUTINE IS BASED ON THE DATAREADER TYPE 9, WHICH WAS TOTALLY
!   CHANGED TO PROVIDE NEW FEATURES AND SKIP TRASH LIKE 'NAMELIST'
!   SOURCECODE USES FORTRAN 90 STANDARD
!   DATA ARE READ INTO DYNAMIC ARRAYS FOR MAXIMUM PERFORMANCE
!   DATAFILES ARE AVAILABLE FOR SEVERAL USERS USING TRNSYS IN A NETWORK
!   DATA TIME SPACING MAY BE > OR < THAN SIMULATION TIME STEP
!   DATA CAN BE INTERPOLATED BY 5_POINT SPLINE
!   TMY, TRY AND VDI 2078 ARE SUPPORTED AS SHORT INPUT OPTION
!
!   THIS ROUTINE IS CALLED BY: LOOPEX
!
!   written by Thomas Lechner, TRANSSOLAR Energietechnik GmbH,
!      March 1996 - February 2001
!************************************************************************************

use TrnsysFunctions ! Use TRNSYS global functions 
use TrnsysConstants ! Use TRNSYS constants
use TrnsysData, ONLY : LUR,LUW,IFORM,LUK ! Use TRNSYS global constants 


use Type109RadiationData
! type109_time_report: 24.02.2005 TW ------------------------------------------------
use type109_time_report

!ak implicit none   ! Force explicit declaration of all variables

! --- TRNSYS declarations -----------------------------------------------------------

!integer :: NI=1, NP=3, ND=0, NO=2, NS=0   ! Nb of inputs, parameters, derivatives, outputs, storage locations

real(8), intent(in)    :: time, xin(*), par(*), t(*)
real(8), intent(inout) :: dtdt(*)
real(8), intent(out)   :: out(*)

integer, intent(inout) :: info(15), iCntrl
real(4)                :: delt,Time0, Tfinal, Tdiff  ! Trnsys timestep and StartTime
! --- END TRNSYS declarations -------------------------------------------------------

! --- Local variables ---------------------------------------------------------------
integer, parameter     :: MAXU=10             !MAXIMUM NUMBER OF UNITS OF TYPE 109 PER SIMULATION
      INTEGER   UNIT(MAXU+1)
      integer :: NPAR, NI, IRADCOL, I_AKI, IC2, I2, IC, JJ, IERR, ILOC
      integer :: ISH, ILATD, ILATM, ILONGD, ILONGM
      integer :: IDP, ILOCB, IDAY, ILOCE, IWARN, NDATA1, NDATA2, DAY 
      integer :: N_LO, NFITLOC, NLASTDAY, NRADI, KBEG, KEND, KKH

character*3, dimension(:), allocatable :: ycheck !CHECK OUTPUT
character*3, dimension(:), allocatable :: ocheck !CHECK OUTPUT
character              :: filename*256
character (len=maxMessageLength) MessErr

integer                :: ifile(5) = 0, updim, readerr
      integer  :: i_low, i_up, n_low(maxu,-1:1)=2 !used to set pointer to actual interpolation

type usave
  character*80         :: stm                 !formatted reading statement
  character*80         :: filename, ofile     !name of opened logical unit
  integer              :: ias                 !data are solar time if ias=1
  integer              :: ifile(5)            !information on data file structure
  integer              :: iform               !formatted reading?
  integer              :: imode               !data reader mode
  integer              :: irad                !number of radiation data columns
  integer              :: ixh                 !repeat data file when exhausted
  integer              :: lu                  !logical unit number for reading
  integer              :: nfit                !size of xfit, yfit
  integer              :: nn                  !number of data to read
  integer              :: nsoll               !number of parameters expected
  integer              :: nupdim              !number of data allocated
  integer              :: sunsetflag          !TIMESCALE UPDATED
  real                 :: ddelt               !time interval of data file
  real                 :: t1                  !time of first data line
  real                 :: tshift              !shift to meet mid point of time interval
  integer, pointer     :: ro(:)               !order for reading of data
  integer, pointer     :: ipol(:)             !interpolation mode
  real,    pointer     :: addf(:)             !addition factor
  real,    pointer     :: samp(:)             !sampling interval before/at/after actual time
  real,    pointer     :: dp(:,:)             !dp contains interpolated values for one day
  real,    pointer     :: prodf(:)            !multiplication factor
  real,    pointer     :: p(:,:)
! p has the following contents:
! (i,-2): integer indicating wether sun is up or down
! (i,-1): time of data point
! (i,0):  same, but sunrise and sunset time replace previous values
! (i,1..): contains values read from data file
  real,    pointer     :: xfit(:)             !workspace for data fitting
  real,    pointer     :: yfit(:)             !workspace for data fitting
  real(8), pointer     :: par(:)              !this compiler forgets include-values
  double precision, pointer :: rres(:)        !result from radiation_processor
end type usave

type (usave), save     :: u(maxu)             !save all data for different calling units
      type (rad_proc_call), save :: r(maxu)              !save all data for different rad procs

      real,    pointer , save    :: xfitloc(:)           !workspace for data fitting
      real,    pointer , save    :: yfitloc(:)           !workspace for data fitting
real                   :: radv(2)
real                   :: eot                 !SEE FUNCTION DEFINITION
real                   :: declin
real                   :: piconv=0.0172028    !2*PI/365.242
real                   :: rdconv=0.0174533    !CONVERT DGREE/RADIANS
real                   :: ecc                 !TIME SHIFT DUE TO ECC. OF EARTH ORBIT
real                   :: sunrise,sunr,sunrmod !TIME FOR ...
real                   :: sunset,suns,sunsmod !TIME FOR ...
real                   :: riseangle           !SUNRISE/SET HOUR ANGLE
real                   :: xtime               !time seems to be intrinsic, so not working with contains
real                   :: shtim(-1:1)         !needed for interpolation
real                   :: xtest               !DUMMY FOR TRIAL READ IN
real                   :: dtk                 !local variable
      real                 :: ttot, timeadd, tanprod, timeloc, tempsum
      real                 :: tn1, tn2, tn2a, tn2e
      real                 :: factkorr, realsum, shloc, xtemp, xlastinterpolation
integer                :: designmonth         !needed for vdi-mode
integer                :: iunit=0             !CURRENT UNIT OF TYPE9 (1 THRU 5)
integer                :: m                   !SAME AS IUNIT, BUT SHORTER
integer                :: ndata               !number of data needed for simulation
integer                :: noutputs            !NUMBER OF OUTPUTS
integer, save          :: falseradflag(maxu)=0  !RADIATION AT NIGHT
integer, save          :: nday                !24 h / delt = # of timesteps per day
integer                :: timep               !POINTER TO TIME ARRAY
integer                :: i,ix, j, k
logical                :: lastdaycut          !SIMULATION ENDS BEFORE END OF DAY
logical                :: frstdaycut          !SIMULATION STARTS LATER THAN DAY
      logical              :: errorflag=.false.    !Errorflag for complete exit
      logical :: true, false

! -------- external procedures ------------------------------------------------------
integer                :: itest,icount,jcount                
real(8)				   :: zwischenzeit(3,10)   
real(8)				   :: summe

!------------------------------------------------------------------------------------
 
!------------------------------------------------------------------------------------
!     Description of Parameters
!         Par(1)  Data Reader Mode /0:user defined normal interpolation
!                 /1: user defined weather data /2: American TMY2 /3: German TRY
!                 /9xx: VDI-2078, dynamische Kuehllastrechnung
!         Par(2)  Number of values (NN) to be read from the data file
!         Par(3)  Data file time interval
!         Par(4)  Time corresponding to first data line
!         Par(5)  rewind if data exhausted 0/1
!         Par(6,7,8,9)Interp. factor, Mult. factor, Add. fact, sampling interv. for first data column
!                     sampling: value before / at / after time   -1 / 0 / 1
!                     Interp. factor/ 0: no / 1: linear / 2: spline / 3: radiation
!          .
!          .
!         Par(2+4*N) Interp. factor for last data column: 0: no / 1: linear / 2: spline
!         Par(3+4*N) Mult factor
!         Par(4+4*N) Add  factor
!         Par(5+4*N) sampling: value before / at / after time   -1 / 0 / 1
!         Par(6+4*N) File logical unit #
!------------------------------------------------------------------------------------
!     Description of Outputs
!         For Mode 0, user defined, no weather data
!          1   user defined variable 1
!          2   user defined variable 2
!          3   ...
!
!         For Mode 1 to 9xx:
!          1   Temp.
!          2   rel.hum.
!          3   wind speed
!          4   wind direction
!		 5   pressure in PA          --  mode 2 (TMY2) / mode 3 (German TRY)
!          5   user defined variable 1 -- only in user defined mode 1
!          6   user defined variable 2 -- only in user defined mode 1
!          7   user defined variable 3 -- only in user defined mode 1
!          8   user defined variable 4 -- only in user defined mode 1
!          9   total extraterrestrial radiation
!          10  ...
!------------------------------------------------------------------------------------

! FUNCTION DEFINITION; EOT USED FOR CALCULATION OF SOLAR TIME ANGLE
eot(x) = -(0.1236*sin(x)   -0.004289*cos(x)+0.1539*sin(2.*x)+0.06078 *cos(2.*x))
ecc(x) = eot(piconv*mod(floor(x/24.),365))

! type109_time_report: 24.02.2005 TW ------------------------------------------------
time_report_unit = 6
time_filter = TFINAL
         time_report = .false.
if (time_report) then
  call get_clock_on_entry('type109')
  zwischenzeit(1,1)=img_high_res_clock(TRUE) ! (start,1 Part)
endif
! ---------------------------------------------------------------------

      if (errorflag) return 1

! --- Initial call to detect the TRNSYS version for which this Type is written -----------------------------------------

if (info(7) .eq. -2) then
  info(12) = 16   ! This component is a TRNSYS 16 Type
  zwischenzeit=0
  itest=0
  return 1
endif

! --- Very last call in simulation --------------------------------------------------
if (info(8).eq.-1) then
! type109_time_report: 24.02.2005 TW ------------------------------
  if (time_report) call get_clock_on_exit('type109')
  if (time_report) call write_clock_time ('type109')
! ---------------------------------------------------------------------
  if (info(8)==-1) then
    if (.not.(header_already_written)) call write_time_report_header
    if (time_report) call write_time_report_summary
    zwischenzeit=zwischenzeit/1000
    Summe=0
    do i=1,10
      write(6,*) 'zwischenzeit',i,' : ',zwischenzeit(3,i)
      Summe=Summe+zwischenzeit(3,i)
    enddo
    write(6,*) 'itest :', itest
    write(6,*) 'Summe : ',Summe
  endif
  return 1    ! Exit - End of the routine for very last call of the simulation
! This block has to be before the normal time steps so we do not do the usual calculations
! outputs were not modified
endif


!!!!!!!!   INITIAL SIMULATION CALL FOR THIS UNIT   !!!!!!!!!
      
xtime = sngl(time) !time is not available in some subroutines, but with assignment it works
! SET UP INFO FOR THIS UNIT OF TYPE
delt  = sngl(getSimulationTimeStep())
time0 = sngl(getSimulationStartTime())
tfinal= sngl(getSimulationStopTime())

if (info(7) == -1) then

  ! SET UP INFO FOR THIS UNIT 
  tdiff =mod(time0,24.)
  if (tdiff>1.) time0=time0-tdiff+1

  iunit = iunit + 1; 
  if (iunit.gt.maxu) then   !TOO MANY TYPE9 UNITS
    call MESSAGES(-1,'The maximum number of instances of Type109 has been exceeded.','fatal',INFO(1),INFO(2))
	return 1
  endif

  unit(iunit) = info(1) ; 
  info(9) = 1; 
  m = iunit
  npar = info(4) ; 
  allocate(u(m)%par(1:npar))
  do i=1,npar; 
    u(m)%par(i)=par(i); 
  enddo

 !check for valid parameter inputs
  select case (int(par(1)))	! check data reader mode (par(1))
    case default
      call TYPECK(-4,INFO,0,1,0)
    case (0)
    case (1:3,911:914,921:924)

    select case (int(par(3)))	! check sky radiation mode (par(3))
      case default
        call TYPECK(-4,INFO,0,3,0)
      case (1:4,99)
    end select
    select case (int(par(4)))	! check  tracking mode (par(4))
      case default
        call TYPECK(-4,INFO,0,4,0)
      case (1:4)
    end select
  end select


  ! TRANSFORM SHORT INPUT OPTION (TMY2,TRY) TO STANDARD INPUT PARAMETER LIST
  u(m)%imode = nint(par(1))
  select case (u(m)%imode)
    case  (0);   u(m)%nn = par(2)            !USER DEFINED READ MODE
    case  (1);   u(m)%nn = 10                !USER DEFINED WEATHER DATA
    case  (2);   u(m)%nn = 7                 !AMERICAN TMY2 MODE
    case  (3);   u(m)%nn = 7                 !GERMAN   TRY  MODE
    case (911:914,921:924)
	  u(m)%nn = 3					!GERMAN   VDI  MODE
    case default; call TYPECK(-4,INFO,0,1,0)
  end select
  if (npar<4) call TYPECK(-4,INFO,0,0,0)			!THERE MUST BE AT LEAST 2 PARAMETERS
  if (u(m)%nn.lt.1.) call TYPECK(-4,INFO,0,0,0)   !NUMBER OF DATA NOT SENSEFULL
  allocate (u(m)%ro(u(m)%nn), u(m)%ipol(1:u(m)%nn), u(m)%prodf(1:u(m)%nn), u(m)%addf(1:u(m)%nn),  &
           u(m)%samp(1:u(m)%nn))
				 u(m)%ro=0.; u(m)%ipol=0.; u(m)%prodf=0.; u(m)%addf=0.; u(m)%samp=0. !cmh 09/05/06 caused error if no wind data was read
  call inip; if (errorflag) return 1        !SET PARAMETERS FROM INPUT
  inquire(unit=u(m)%lu,name=u(m)%filename)

  hdelt = delt/2.; cdelt = delt*0.01         !USED BY RADIATION PROCESSOR
  r(m)%infloc(1)=info(1); r(m)%infloc(2)=info(2)        !USED BY RADIATION PROCESSOR


  ! LU MAY NOT BE LESS THAN 10
  if (u(m)%lu.lt.10) then
    CALL MESSAGES(215,'','fatal',INFO(1),INFO(2))   
    return 1
  endif

  ! DATA FILE TESTING
  call infiletest; if (errorflag) return 1
  call input_init; if (errorflag) return 1 !INPUTS IN CASE OF RADIATION PROCESSING

  ni = info(3)
  info(6) = noutputs                !PICK UP NUMBER OF OUTPUTS
  call typeck(1,INFO,NI,-1,0)       !DON'T CHECK NUMBER OF PARAMETERS
                                            !AND SET NUMBER OF OUTPUTS
  if (u(m)%nsoll/=npar) call TYPECK(-4,INFO,0,0,0)    !CHECK NUMBER OF PARAMETERS

  allocate (ycheck(info(3)));allocate (ocheck(info(6)))
  if (u(m)%imode==0) then
    ocheck='DM1'
  else
    ocheck(1:4)=(/'TE1','PC1','VE1','DG1'/); ocheck(5:8)='DM1'
    ocheck(9:11)=(/'IR1','DG1','DG1'/)
    do i=2,r(m)%nsurf+2
      ocheck(6*i:6*i+5)=(/'IR1','IR1','IR1','IR1','DG1','DG1'/)
    enddo
  endif

  call rcheck(info,ycheck,ocheck)


  ! ALLOCATE ARRAYS AND SET POINTER TO CORRESPONDING IUNIT
  nday = nint(24./delt)
  u(m)%nfit = nint(24./u(m)%ddelt)
  if (u(m)%ixh/=0) then                           !REWIND OPTION ACTIVE
    ifile(3)=ifile(2)+ifile(3); ifile(2)=0; ifile(5)=0
  endif

  updim=ifile(3)+4; u(m)%nupdim = updim			!all data are read into u%p here
  allocate (u(m)%dp(nday,u(m)%irad))
  allocate (u(m)%p(updim,-2:u(m)%nn)); 	
  do icount = 1,updim
    do jcount =-2,u(m)%nn
      u(m)%p(icount,jcount)=0	
    enddo
  enddo

  allocate (u(m)%xfit(u(m)%nfit), u(m)%yfit(u(m)%nfit))
  allocate (xfitloc(u(m)%nfit), yfitloc(u(m)%nfit))

  ! READ THE DATA FROM FILE TO ARRAY, read and close data file
  call dataread; if (errorflag) return 1; close(u(m)%lu)

!initialize outputs
  out(1:noutputs) = 0.d0
!initialize sunrise and sunset angles
  call update_time_sun_radiation

  return 1

endif
!!!!!!!   END OF INITIAL SIMULATION CALL FOR THIS UNIT   !!!!!!!!!
!-----------------------------------------------------------------

! ITERATIVE CALL -------------------------------------------------

      do iunit=1,maxu+1        ! FIGURE OUT WHICH OF THE TYPE9 UNITS THIS IS
        if(unit(iunit) == info(1)) exit
      enddo
      if (iunit>maxu)  then
	    call MESSAGES(216,'','fatal',INFO(1),INFO(2))
		return 1
	  endif
      m = iunit


! Read the slope and azimuth
	  if (PAR(4) ==  1) then
        do i=1,r(m)%nsurf
          r(m)%slopes(i) = xin(2*i)
          r(m)%axes(i)  = xin(2*i+1)
        enddo
	  endif

! UPDATE TIME, SUNRISE, SUNSET, RADIATION
      call update_time_sun_radiation

! type109_time_report: 24.02.2005 TW ------------------------------
       if (time_report) then
        zwischenzeit(2,1)=img_high_res_clock(FALSE) ! (End,1 Part)
          if (time>tfinal/10) Then
	 	     zwischenzeit(1,2)=img_high_res_clock(TRUE) ! (Start,2 Part)
          else
		     zwischenzeit(1,10)=img_high_res_clock(TRUE) ! (Start,10 Part)
          endif
        zwischenzeit(3,1)=zwischenzeit(3,1)+zwischenzeit(2,1)-zwischenzeit(1,1) ! (dT,1 Part)
       endif
! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
 
! INTERPOLATE DATA
!WHICH ARRAY ELEMENT CONTAINS ACTUAL DATA FOR SAMPLING -1/0/1 ?
      do i=-1,1,1
        shtim(i)=xtime+(u(m)%tshift-i+1)*u(m)%ddelt*0.5

        do while (u(m)%p(n_low(m,i)+1,-1)<shtim(i))
         n_low(m,i)=n_low(m,i)+1
		 itest =itest+1
	    enddo

        if (abs(u(m)%p(n_low(m,i)-1,-1)-shtim(i))<abs(u(m)%p(n_low(m,i),-1)-shtim(i))) n_low(m,i)=n_low(m,i)-1
      enddo
      iradcol=0

! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
! ---------------------------------------------------------------------
! type109_time_report: 24.02.2005 TW ------------------------------
       if (time_report) then
          if (time>tfinal/10) Then
        zwischenzeit(2,2)=img_high_res_clock(FALSE) ! (End,2 Part)
        zwischenzeit(3,2)=zwischenzeit(3,2)+zwischenzeit(2,2)-zwischenzeit(1,2) ! (dT,2 Part)
          else
        zwischenzeit(2,10)=img_high_res_clock(FALSE) ! (End,10 Part)
        zwischenzeit(3,10)=zwischenzeit(3,10)+zwischenzeit(2,10)-zwischenzeit(1,10) ! (dT,10 Part)
          endif
        zwischenzeit(1,3)=img_high_res_clock(TRUE) ! (Start,3 Part)
       endif
! ---------------------------------------------------------------------

      do i=1,u(m)%nn
        k = u(m)%samp(i); dtk = (k-u(m)%tshift)*u(m)%ddelt*0.5
        select case (u(m)%ipol(i))

         case (0)                            !NO INTERPOLATION
           if (k==0) then
             if (xtime-u(m)%p(n_low(m,k),0)-dtk > u(m)%ddelt/2.) then
               out(i)=u(m)%p(n_low(m,k)+1,i)
             else
               out(i)=u(m)%p(n_low(m,k),i)
             end if
           else
             out(i)=u(m)%p(n_low(m,k),i)
           end if

         case (1)                            !LINEAR INTERP.
            if (u(m)%p(n_low(m,k),-1)>=xtime-dtk) n_low(m,k) = n_low(m,k)-1
           i_low = n_low(m,k); i_up  = i_low+1; i_aki = 2
           out(i)=u(m)%p(i_up,i)+(xtime-dtk-delt/2.-u(m)%p(i_up,-1))/u(m)%ddelt *(u(m)%p(i_up,i)-u(m)%p(i_low,i))

         case (2)                            !NATURAL SPLINE
           i_low = max(1,n_low(m,k)-3)
           i_up  = min(n_low(m,k)+3,updim)
           i_aki = i_up-i_low+1
           ic2 = 1
           do i2=i_low,i_up
             u(m)%xfit(ic2) = u(m)%p(i2,-1)
             u(m)%yfit(ic2) = u(m)%p(i2,i)
             ic2=ic2+1
           enddo
           call akimaip(luw,i_aki,u(m)%xfit,u(m)%yfit,1,xtime-delt/2.-dtk,radv)

           out(i) = radv(1)
! type109_time_report: 24.02.2005 TW ------------------------------
       if (time_report) then
        zwischenzeit(2,3)=img_high_res_clock(FALSE) ! (End,3 Part)
        zwischenzeit(1,4)=img_high_res_clock(TRUE) ! (Start,4 Part)
        zwischenzeit(3,3)=zwischenzeit(3,3)+zwischenzeit(2,3)-zwischenzeit(1,3) ! (dT,3 Part)
       endif
! ---------------------------------------------------------------------

         case (3) !RADIATION INTERPOLATION TREATING DATA AS INTEGRATED OVER TIMESTEP
           iradcol = iradcol+1
           if (iradcol>u(m)%irad) cycle; ic = nint(mod(xtime,24.)/delt)
           if (ic==0) ic = nint(24./delt)
           out(i) = u(m)%dp(ic,iradcol)

        end select
      enddo
! --- let t_amb and rh: 25.10.2005 YY ---------------------------------
      r(m).t_amb = out(1)
	  r(m).rh_amb = out(2)

           !CALL RADIATION PROCESSOR
      select case (u(m)%imode)
        case (1)                              !USER DEFINED WEATHER DATA
          select case (r(m)%radmode)
            case (4,5)						  !I and T/Hum
              r(m)%iges   = out(5)
            case (6)						  !Ib and Id
              r(m)%ib     = out(5); r(m)%id   = out(6)
            case (7)						  !I adn Ibn
              r(m)%iges   = out(5); r(m)%ibdn = out(6)
            case (8)						  !I and Id
              r(m)%iges   = out(5); r(m)%id   = out(6)
            case (9)						  !Ibn and Id
              r(m)%ibdn     = out(5); r(m)%id   = out(6)
          end select
          out(5:6)=0.
          do jj=        5,u(m)%nn-2; out(jj)=out(jj+2); end do
          do jj=u(m)%nn-1,u(m)%nn;   out(jj)=0.;        end do
        case (2)                              !TMY2-MODE
          r(m)%iges = out(5); r(m)%id= out(6); out(5)=out(7); out(6:8)=0.
        case (3)                              !TRY-MODE
          r(m)%ib   = out(5); r(m)%id= out(6); out(5)=out(7); out(6:8)=0.
        case (9:)                             !VDI-MODE
          r(m)%iges = out(2); r(m)%id= out(3); out(2:8)=0.
      end select

      if (r(m)%skymode==99) then          !zero radiation
	  r(m)%iges=0.;r(m)%ib=0.;r(m)%id=0.
	  endif

      call rad_proc(xtime, r(m), u(m)%rres)
! type109_time_report: 24.02.2005 TW ------------------------------
       if (time_report) then
        zwischenzeit(2,4)=img_high_res_clock(FALSE) ! (End,4 Part)
        zwischenzeit(1,5)=img_high_res_clock(TRUE) ! (Start,5 Part)
        zwischenzeit(3,4)=zwischenzeit(3,4)+zwischenzeit(2,4)-zwischenzeit(1,4) ! (dT,4 Part)
       endif
! ---------------------------------------------------------------------

      if (u(m)%imode>0) out(9:r(m)%upper+8) = u(m)%rres(1:r(m)%upper)

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('type109')
      if (time_report) call write_clock_time ('type109')
! ---------------------------------------------------------------------

      return 1


      contains    !following subroutines until 'END SUBROUTINE TYPE109' are
                  ! only locally available in subroutine TYPE109 and share
                  ! all values
!*****************************************************************************
      subroutine inip

! initialize parameters for data-reading
      integer   ntmy2, ntry, nvdi
      parameter (ntmy2=33, ntry=33, nvdi=20)
      real   :: tmy2(ntmy2), try(ntry), vdi(nvdi)         !initialization values

      type (readvar)          :: var(50), v
      character               :: instr*200, part*8, dirw13*13, solar*5, LineStr
	  character(len=200)	  :: ErrMsg1
      integer                 :: j, pos, lcount, c, indx, indx2, cusrv
      integer                 :: maxcol, mincol, mincount, fcount
      integer                 :: nskip(50), nu(50), iounit
      real                    :: long, dum(100), vdum(100), ttest
      logical                 :: fileopen=.TRUE.
      logical                 :: firstrad


      data tmy2 /2, 8, 1, 1, 1,  2,0.1,0,0,  1,1,0,0,  1,0.1,0,0,  1,1,0,0,  3,3.6,0,-1,  3,3.6,0,-1,  1,100,0,0/
! TMY2:    TEMP    RHUM      WIND    WDIR I         Id     press
! every variable needs the four values: interp, multipl. add. sampling
      data try  /3, 8, 1, 1, 1, 2, 1.,0,0,1,100,0,0,1,1,0,0,1,1,0,0,3,3.6,0,-1,3,3.6,0,-1,1,100,0,0/
!  TRY:    TEMP      RHUM    WIND    WDIR         Ib         Id     press
      data vdi /0,9,1,0,1, 1,1,0,0, 3,3.6,0,0, 3,3.6,0,0, 50,0,1/
 
! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('inip')
      endif
! ---------------------------------------------------------------------

      u(m)%irad=0; readerr=0; part=' '; lcount=0; c=0; maxcol=0
      indx=0; indx2=0; nu=0; cusrv=0; fcount=6
      select case (u(m)%imode)
       case (0)                               !NORMAL USER DEFINED INPUT
         u(m)%iform = 0                       !NO RADIATION
         u(m)%tshift= 0.
         do i=1,u(m)%nn; u(m)%ro(i)=i; enddo
         u(m)%ddelt=u(m)%par(3); u(m)%t1=u(m)%par(4)
         u(m)%irad  = 0; u(m)%ixh=u(m)%par(5)
         do i=1,u(m)%nn; u(m)%ipol(i) = u(m)%par(4*i+2)
           if (u(m)%ipol(i)<0.or.u(m)%ipol(i)>2) call TYPECK(-4,INFO,0,4*i+2,0)
           u(m)%prodf(i) = u(m)%par(4*i+3)
           u(m)%addf(i)  = u(m)%par(4*i+4)
           u(m)%samp(i)  = u(m)%par(4*i+5)
         enddo
         u(m)%lu    = nint(u(m)%par(4*u(m)%nn+6))
         inquire(unit=u(m)%lu,name=u(m)%filename)
         u(m)%nsoll = 4*u(m)%nn+6

       case (1)                               !USER DEFINED INPUT WITH RADIATION PROCESSOR CALL
         u(m)%iform       = 0                 !free format reading
         u(m)%lu          = u(m)%par(2)
         inquire(unit=u(m)%lu,name=u(m)%filename)
         r(m)%skymode     = u(m)%par(3)           !skymodel
         r(m)%trackmode   = u(m)%par(4)           !trackingmode
         u(m)%ixh         = 1
         u(m)%nsoll       = 4
         firstrad=.FALSE.

      !data file check
         read(u(m)%lu,'(A)',err=771) dirw13
         if (dirw13.ne.'<userdefined>') then
		   call MESSAGES(217,'','fatal',INFO(1),INFO(2))
		   return
		 endif

         do while (part /= '<data>')                          !read data file header
           read(u(m)%lu,'(A)')instr
           pos = index(instr,">")
           read(instr(1:pos),*) part
           lcount=lcount+1
           instr =instr(pos+1:)       !cut first keyword <abc>

           select case (part)     !read definitions
            case ('<longitu')
              read(instr,*,err=334) long
              long=-long
            case ('<latitud')
              read(instr,*,err=334) r(m)%latit
            case ('<interva')
              read(instr,*,err=334) u(m)%ddelt
            case ('<gmt>')
              instr = adjustl(instr)
              read(instr,'(A5)',err=334) solar
              read(instr,*,err=310) r(m)%gmt
              r(m)%shift=-r(m)%gmt*15.+long
              cycle
  310         if (solar == "solar") then  !if gmt is "solar" instead of a number, treat user def data as solar time
                r(m)%ias = 1
              else
                goto 334
              endif
            case ('<firstti')
              read(instr,*,err=334) u(m)%t1
            case ('<var>')
              read(instr,*,err=334) v%name,part,v%icol,part,v%ipol,part,v%add,part,v%mult,part,v%samp
              if (v%icol<1) cycle
              select case (v%name)
                case ('TAMB','RHUM','WSPEED','WDIR')
                  v%pname='==> '//v%name
                case ('IBEAM_H','IGLOB_H','IBEAM_N','IDIFF_H')
                  if (.not.firstrad.and.v%icol/=0) then
                    firstrad=.TRUE.; u(m)%tshift = v%samp
                  else
                    v%samp = u(m)%tshift
                  endif
                  v%pname='==> '//v%name
                case default
                  v%pname='    '//v%name
                  if (v%icol>0) cusrv = cusrv+1
              end select
              c=c+1
              if (cusrv>4) goto 338       !not more than 4 usrdefined data columns
              var(c)=v
            case ('<data>')
              goto 330                !end of header section
            case default
              goto 334                !found wrong keyword
           end select

           if (lcount>100) goto 333
         end do
         return

!available radiation data analysis
330     do j=1,c
          if(var(j)%name=='IBEAM_N')indx=indx+1
          if(var(j)%name=='IBEAM_H')indx=indx+2
          if(var(j)%name=='IDIFF_H')indx=indx+4
          if(var(j)%name=='IGLOB_H')indx=indx+8
        end do
        select case (indx)
          case (1:4); goto 335    !Not enough data

		  case (5)				  !Ibn and Id are inputs
				r(m)%radmode=9
          case (6)                !Ib and Id  are inputs
                r(m)%radmode=6
          case (7); goto 336      !Modify data
          case (8)                !Iglob: reduced Reindl or full Reindl
                do j=1,c
                  if(var(j)%name=='TAMB')indx2=indx2+1
                  if(var(j)%name=='RHUM')indx2=indx2+2
                end do
         if (indx2==3) then; r(m)%radmode=5; else; r(m)%radmode=4;endif
          case (9)                !I  and Ibn are inputs
                r(m)%radmode=7
          case (10); goto 337     !Modify data
          case (11)               !I, Ib  and Ibn  are inputs; kick off Ib
                r(m)%radmode=7
                doloop11:do j=1,c
                  if (var(j)%name=='IBEAM_H') then
                    do k=j,c-1
                      var(k) = var(k+1)
                    end do
                    c=c-1; exit doloop11
                  endif
                end do doloop11
          case (12)               !I  and Id  are inputs
                r(m)%radmode=8
          case (13)               !I, Id  and Ibn  are inputs; kick off Id
                r(m)%radmode=7
                doloop13:do j=1,c
                  if (var(j)%name=='IDIFF_H') then
                    do k=j,c-1
                      var(k) = var(k+1)
                    end do
                    c=c-1; exit doloop13
                  endif
                end do doloop13
          case (14)               !I, Id  and Ib  are inputs; kick off Ib
                r(m)%radmode=8
                doloop14:do j=1,c
                  if (var(j)%name=='IBEAM_H') then
                    do k=j,c-1
                      var(k) = var(k+1)
                    end do
                    c=c-1; exit doloop14
                  endif
                end do doloop14
          case (15)               !I, Id, Ib  and Ibn  are inputs; kick off Ib and Id
                r(m)%radmode=7
                doloop151:do j=1,c
                  if (var(j)%name=='IBEAM_H') then
                    do k=j,c-1
                      var(k) = var(k+1)
                    end do
                    c=c-1; exit doloop151
                  endif
                end do doloop151
                doloop152:do j=1,c
                  if (var(j)%name=='IDIFF_H') then
                    do k=j,c-1
                      var(k) = var(k+1)
                    end do
                    c=c-1; exit doloop152
                  endif
                end do doloop152
          case default
                goto 334
        end select

!print out how data are interpreted
        call MESSAGES(-1,'','Notice',INFO(1),INFO(2))    
        write(luw,*)' VARIABLE              COLUMN # IPOL     MULTIPLIER    ADD_FACTOR   SAMPLING'
        write(luw,*)'----------------------------------------------------------------------------'
        do j=1,c
          select case (var(j)%name)
            case ('IBEAM_H','IGLOB_H','IBEAM_N','IDIFF_H')
              write(luw,'(2X,A20,3X,I2,7X,A2,7X,F9.4,5X,F9.4,8X,I2)') var(j)%pname,var(j)%icol,' -',var(j)%mult,var(j)%add,var(j)%samp
            case default
              write(luw,'(2X,A20,3X,I2,7X,I2,7X,F9.4,5X,F9.4,8X,I2)') var(j)%pname,var(j)%icol,var(j)%ipol,var(j)%mult,var(j)%add,var(j)%samp
          end select
        end do

        write(luw,*)'----------------------------------------------------------------------------'
        write(luw,'(A,F7.2)')'  LONGITUDE               ',-long
        write(luw,'(A,F7.2)')'  LATITUDE                ',r(m)%latit
        write(luw,'(A,I3)')  '  TIME SHIFT TO GMT        ',r(m)%gmt
        write(luw,'(A,F7.2)')'  DATA TIME STEP          ',u(m)%ddelt
        write(luw,'(A,F7.2)')'  STARTTIME FOR DATA      ',u(m)%t1
        write(luw,'(A,F7.2)')'  SAMPLING FOR RAD. DATA  ',u(m)%tshift
        write(luw,*)'----------------------------------------------------------------------------'
        select case (r(m)%radmode)
          case(4)
            write(luw,*)' RADIATION MODE: REDUCED REINDL'
            u(m)%irad = 1
          case(5)
            write(luw,*)' RADIATION MODE: FULL REINDL'
            u(m)%irad = 1
          case(6)
            write(luw,*)' RADIATION MODE: IBEAM_H AND IDIFF_H ARE INPUTS'
            u(m)%irad = 2
          case(7)
            write(luw,*)' RADIATION MODE: IGLOB_H AND IBEAM_N ARE INPUTS'
            u(m)%irad = 2
          case(8)
            write(luw,*)' RADIATION MODE: IGLOB_H AND IDIFF_H ARE INPUTS'
            u(m)%irad = 2
          case(9)
            write(luw,*)' RADIATION MODE: IBEAM_N AND IDIFF_H ARE INPUTS'
            u(m)%irad = 2

        end select
        write(luw,*)

!new reading order due to variable name
        u(m)%samp(1:u(m)%nn) = 0.d0
        do j=1,c
          select case (var(j)%name)
            case('TAMB')
                  nu(1)=var(j)%icol
                  u(m)%prodf(1) = var(j)%mult
                  u(m)%addf(1)  = var(j)%add
                  u(m)%samp(1)  = var(j)%samp
                  u(m)%ipol(1)  = var(j)%ipol
            case('RHUM')
                  nu(2)=var(j)%icol
                  u(m)%prodf(2) = var(j)%mult
                  u(m)%addf(2)  = var(j)%add
                  u(m)%samp(2)  = var(j)%samp
                  u(m)%ipol(2)  = var(j)%ipol
            case('WSPEED')
                  nu(3)=var(j)%icol
                  u(m)%prodf(3) = var(j)%mult
                  u(m)%addf(3)  = var(j)%add
                  u(m)%samp(3)  = var(j)%samp
                  u(m)%ipol(3)  = var(j)%ipol
            case('WDIR')
                  nu(4)=var(j)%icol
                  u(m)%prodf(4) = var(j)%mult
                  u(m)%addf(4)  = var(j)%add
                  u(m)%samp(4)  = var(j)%samp
                  u(m)%ipol(4)  = var(j)%ipol
            case('IBEAM_H','IGLOB_H')
				  nu(5)=var(j)%icol
                  u(m)%prodf(5) = var(j)%mult*3.6
                  u(m)%addf(5)  = var(j)%add
                  u(m)%samp(5)  = var(j)%samp
                  u(m)%ipol(5)  = var(j)%ipol
            case('IBEAM_N','IDIFF_H')
				  if ( (r(m)%radmode == 9) .and. (var(j)%name == 'IBEAM_N'  ) ) then 
					!Case when Ibn and Id are inputs
				    nu(5)=var(j)%icol
                    u(m)%prodf(5) = var(j)%mult*3.6
                    u(m)%addf(5)  = var(j)%add
                    u(m)%samp(5)  = var(j)%samp
                    u(m)%ipol(5)  = var(j)%ipol
				  else
					!Other cases
                  nu(6)=var(j)%icol
                  u(m)%prodf(6) = var(j)%mult*3.6
                  u(m)%addf(6)  = var(j)%add
                  u(m)%samp(6)  = var(j)%samp
                  u(m)%ipol(6)  = var(j)%ipol
				  endif
            case default; fcount=fcount+1; nu(fcount)=var(j)%icol
                  u(m)%prodf(fcount) = var(j)%mult
                  u(m)%addf(fcount)  = var(j)%add
                  u(m)%samp(fcount)  = var(j)%samp
                  u(m)%ipol(fcount)  = var(j)%ipol
          end select
        end do
        do j=1,10; if (nu(j)==0) nu(j)= 11; enddo


!prepare free format reading
        mincol = var(1)%icol; maxcol = var(1)%icol
        do j=1,c-1            !sort <var> to ascending column #
          mincol   = var(j)%icol
          mincount = j
          maxcol = max(maxcol,var(j+1)%icol)
          do k=j,c-1
            if (mincol>var(k+1)%icol) then
              mincol   = var(k+1)%icol
              mincount = k+1
            end if
          end do
          var(11) = var(j)
          var(j)  = var(mincount)
          var(mincount)  = var(11)
        end do

!open temporary internal file
        iounit=899
        do while (fileopen)
          inquire(iounit, opened=fileopen)
          iounit = iounit+1
        end do
        open(iounit, status='scratch')

        vdum(11)=0.
  	    u(m)%ifile(1) = lcount+1
        do                    !read desired data, skip the rest
          lcount=lcount+1
          read(u(m)%lu,*,end=340,err=780)(vdum(j),j=1,maxcol)
          write(iounit,'(10(e12.6,1X))')(vdum(nu(j)),j=1,10)
        end do

340     close(u(m)%lu); u(m)%lu = iounit  !change lunit #
        u(m)%nn=fcount
        u(m)%ipol(5) = 3; u(m)%ipol(6) = 3 !set radiation to interp mode 3
        rewind(u(m)%lu)
        do i=1, u(m)%nn; u(m)%ro(i)=i; end do

       case (2)                                 !TMY2-PARAMETERS
             ! DATA ORDER: TEMP RHUM WIND WDIR ITOT IDIF PRES 1 2 3 4 5 6 7
             ! READ ORDER: ITOT IDIF TEMP RHUM PRES WDIR WIND 5 6 1 2 7 4 3
         u(m)%ro(1) =5; u(m)%ro(2) =6; u(m)%ro(3)=1; u(m)%ro(4) =2
         u(m)%ro(5) =7; u(m)%ro(6) =4; u(m)%ro(7) =3
         u(m)%stm = '(17X,F4.0,8X,F4.0,34X,F4.0,8X,F3.0,2X,F4.0,2X,F3.0,2X,F3.0)'
         u(m)%iform=1
         u(m)%ddelt=tmy2(3); u(m)%t1=tmy2(4); u(m)%tshift=-1
         u(m)%ixh=tmy2(5)
         do i=1,u(m)%nn
	         u(m)%ipol(i) = tmy2(4*i+2)
              u(m)%prodf(i) = tmy2(4*i+3)
               u(m)%addf(i) = tmy2(4*i+4)
               u(m)%samp(i) = tmy2(4*i+5)
         enddo
         u(m)%lu    = u(m)%par(2)
         inquire(unit=u(m)%lu,name=u(m)%filename)
         u(m)%irad  = 2
         u(m)%nsoll = 4
         u(m)%ixh   = 1
         r(m)%radmode   = 8      !radiationmode I_global_normal, I_diffuse_normal are inputs
         r(m)%skymode   = u(m)%par(3) !skymodel
         r(m)%trackmode = u(m)%par(4) !trackingmode

       case (3)                                 !TRY-PARAMETERS
             ! DATA ORDER: TEMP RHUM WIND WDIR IBEAM IDIF PRES  1 2 3 4 5 6 7
             ! READ ORDER: WDIR WIND PRES TEMP RHUM IBEAM IDIF  4 3 7 1 2 5 6
         u(m)%ro(1) =4; u(m)%ro(2) =3; u(m)%ro(3)=7; u(m)%ro(4) =1
         u(m)%ro(5) =2; u(m)%ro(6) =5; u(m)%ro(7) =6
         u(m)%stm = '(13X,F3.0,X,F4.1,14X,F6.1,1X,F5.1,X,F4.2,X,E4.0,X,E4.0)'
         u(m)%iform=1
         u(m)%ddelt=try(3); u(m)%t1=try(4); u(m)%tshift=-1
         u(m)%ixh=try(5)
         do i=1,u(m)%nn
	          u(m)%ipol(i) = try(4*i+2)
               u(m)%prodf(i) = try(4*i+3)
                u(m)%addf(i) = try(4*i+4)
                u(m)%samp(i) = try(4*i+5)
         enddo
         u(m)%lu    = u(m)%par(2)
         inquire(unit=u(m)%lu,name=u(m)%filename)
         u(m)%irad  = 2
         u(m)%nsoll = 4
         u(m)%ixh   = 1
         r(m)%radmode   = 6      !radiationmode I_beam_normal, I_diffuse_normal are inputs
         r(m)%skymode   = u(m)%par(3) !skymodel
         r(m)%trackmode = u(m)%par(4) !trackingmode

       case (9:)                                 !VDI-PARAMETERS
             ! DATA ORDER: IDIR IDIF TEMP  1 2 3 0 0 0
             ! READ ORDER: TEMP IDIR IDIF  3 1 2
         u(m)%ro(1) =2; u(m)%ro(2) =3; u(m)%ro(3)=1

         ttest = mod(time0,24.);
         if (ttest.ne.0.) goto 339
         u(m)%iform=1
         iloc = u(m)%imode-900
         if (iloc<20) then                        ! 91x means truebes Wetter
              u(m)%stm = '(10X,F10.0,F10.0,20X,'
              iloc = iloc -10
           else                                   ! 92x means heiteres Wetter
              u(m)%stm = '(30X,F10.0,F10.0,'
              iloc = iloc -20
         end if

         select case (iloc)
           case (1)                               ! 9x1 means Kuehllastzone 1
              u(m)%stm = trim(u(m)%stm)//'F10.1)'
           case (2)                               ! 9x2 means Kuehllastzone 2
              u(m)%stm = trim(u(m)%stm)//'10X,F10.1)'
           case (3)                               ! 9x3 means Kuehllastzone 3
              u(m)%stm = trim(u(m)%stm)//'20X,F10.1)'
           case (4)                               ! 9x4 means Kuehllastzone 4
              u(m)%stm = trim(u(m)%stm)//'30X,F10.1)'
         end select

         u(m)%ddelt=vdi(3)
         select case (ceiling(time0/24))  !VDI needs a     day repeatedly
           case (:31);        r(m)%designday =  24; designmonth =  1
           case (32:59);      r(m)%designday =  51; designmonth =  2
           case (60:90);      r(m)%designday =  81; designmonth =  3
           case (91:120);     r(m)%designday = 110; designmonth =  4
           case (121:151);    r(m)%designday = 141; designmonth =  5
           case (152:181);    r(m)%designday = 172; designmonth =  6
           case (182:212);    r(m)%designday = 204; designmonth =  7
           case (213:243);    r(m)%designday = 236; designmonth =  8
           case (244:273);    r(m)%designday = 260; designmonth =  9
           case (274:304);    r(m)%designday = 296; designmonth = 10
           case (305:334);    r(m)%designday = 324; designmonth = 11
           case (335:365);    r(m)%designday = 356; designmonth = 12
         end select
         u(m)%t1=floor(time0/24)*24
         u(m)%tshift=0; u(m)%ixh=vdi(5)
         do i=1,u(m)%nn
               u(m)%ipol(i) = vdi(4*i+2)
              u(m)%prodf(i) = vdi(4*i+3)
               u(m)%addf(i) = vdi(4*i+4)
               u(m)%samp(i) = vdi(4*i+5)
         enddo
         u(m)%lu        = u(m)%par(2)
         inquire(unit=u(m)%lu,name=u(m)%filename)
         u(m)%irad      = 4
         u(m)%nsoll     = 4
         r(m)%ivdi      = 1
         r(m)%latit     = vdi(18)
         r(m)%shift     = vdi(19)
         u(m)%ias       = vdi(20)
         r(m)%radmode   = 8      !radiationmode I_global_normal, I_diffuse_normal are inputs
         r(m)%skymode   = u(m)%par(3) !skymodel
         r(m)%trackmode = u(m)%par(4) !trackingmode

       case default
        write(luw,*)'Choice of Mode in TYPE 109 (Datareader) not valid'
	  stop
      end select


      if (u(m)%ias==1) then
         r(m)%shift = 0.    !TREAT AS SOLAR
         r(m)%ias    = 1
      endif
      u(m)%ofile = u(m)%filename


      return

333   call MESSAGES(218,'','fatal',INFO(1),INFO(2)); return
334   call MESSAGES(219,'','fatal',INFO(1),INFO(2)); return
335   call MESSAGES(220,'','fatal',INFO(1),INFO(2)); return
336   call MESSAGES(221,'','fatal',INFO(1),INFO(2)); return
337   call MESSAGES(222,'','fatal',INFO(1),INFO(2)); return
338   call MESSAGES(223,'','fatal',INFO(1),INFO(2)); return
339   call MESSAGES(224,'','fatal',INFO(1),INFO(2)); return
771   call MESSAGES(217,'','fatal',INFO(1),INFO(2)); return

780   WRITE(LineStr,*) lcount + 1
      ErrMsg1 = 'Error ocurred in line:'//TRIM(ADJUSTL(LineStr))//'.'
      call MESSAGES(214,ErrMsg1,'fatal',IUNIT,ITYPE)   
	  return

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('inip')
      if (time_report) call write_clock_time ('inip')
! ---------------------------------------------------------------------
      end subroutine inip
!*****************************************************************************



!*****************************************************************************
      subroutine input_init
!*****************************************************************************
! DETERMINE NUMBER OF INPUTS (GROUND REFLECTANCE, SLOPE, AZIMUTH)
! CHECK FOR ERRORS

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('input_init')
      endif
! ---------------------------------------------------------------------

      r(m)%ninp  = info(3)                                            !# OF INPUTS
      r(m)%nsurf = (r(m)%ninp-1)/2;r(m)%nsurf = max(r(m)%nsurf,1) !# OF SURFACES
      allocate (r(m)%slopes(r(m)%nsurf))  !SHAPE ARRAY CORRESPONDING TO # OF SURFACES
      allocate (r(m)%axes(r(m)%nsurf))    !SHAPE ARRAY CORRESPONDING TO # OF SURFACES
      r(m)%upper = 6*r(m)%nsurf+9         !NUMBER OF OUTPUTS FROM RAD_PROC
      allocate (u(m)%rres(r(m)%upper))    !SHAPE ARRAY CORRESPONDING TO # OF SURFACES
      r(m)%rho= xin(1)
      do i=1,r(m)%nsurf
          r(m)%slopes(i) = xin(2*i)
          r(m)%axes(i)  = xin(2*i+1)
      enddo
      r(m)%sinlat=sin(r(m)%latit)
      r(m)%coslat=cos(r(m)%latit)
      r(m)%tanlat=tan(r(m)%latit)

! DETERMINE NUMBER OF OUTPUTS
      if (u(m)%imode==0) noutputs = u(m)%nn
      if (u(m)%imode>0)  noutputs = 17+6*r(m)%nsurf

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('input_init')
      if (time_report) call write_clock_time ('input_init')
! ---------------------------------------------------------------------

      end subroutine input_init
!*****************************************************************************


!******************************************************************************
      subroutine infiletest
!******************************************************************************
! CHECKS WETHER DATAFILE IS OK AND COUNTS NUMBER OF LEADING COMMENT LINES
! AND NUMBER OF DATA LINES.
! WHAT HAPPENS HERE?
!  readerr=-1    DATA File empty
!  readerr=1     DATA LINE FOLLOWED BY COMMENT LINE; STOP SIMULATION
!  readerr=2     MORE THAN maxcom COMMENT LINES; STOP SIMULATION
!  readerr=3     u(m)%ixh = 0, BUT DATA EXHAUSTED
!  readerr=4     time0 < t1
!  readerr=5     number of data lines does not fit a multiple of one day for
!                repeated rewind



!  ifile(1)      number of comment lines; 20 lines maximum
!  ifile(2)      number of skipped lines from first data line to line 
!                corresponding to time0
!  ifile(3)      number of data lines which could be read in
!  ifile(4)      number of data lines which should be read in
!  ifile(5)      more data lines


      real       :: x
      integer    :: maxcom=20
      integer    :: nskip, xflag=0, cc=0, ix
      character  :: dirbuch*1, dirword*5, cdum*80, dirLat*1, LineStr
	  character(len=200)	  :: ErrMsg1

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('infiletest')
      endif
! ---------------------------------------------------------------------

      ndata=0; nskip=0; readerr=0; xflag=0; cc=0
      nskip = max(ceiling((time0-u(m)%t1)/u(m)%ddelt),0)
      ndata  = nint((tfinal-time0)/u(m)%ddelt)
      ifile=0; ifile(4) = ndata
	
      goto 201    !Jump over ...


770   call MESSAGES(225,'','fatal',INFO(1),INFO(2)); return
772   call MESSAGES(226,'','fatal',INFO(1),INFO(2)); return
773   call MESSAGES(227,'','fatal',INFO(1),INFO(2)); return
779   call MESSAGES(228,'','fatal',INFO(1),INFO(2)); return
801   call MESSAGES(229,'','fatal',INFO(1),INFO(2)); return
802   call MESSAGES(230,'','fatal',INFO(1),INFO(2)); return
803   call MESSAGES(231,'','fatal',INFO(1),INFO(2)); return
804   call MESSAGES(232,'','fatal',INFO(1),INFO(2)); return
805   call MESSAGES(233,'','fatal',INFO(1),INFO(2)); return
806   call MESSAGES(234,'','fatal',INFO(1),INFO(2)); return


780   WRITE(LineStr,*) cc
      ErrMsg1 = 'Error ocurred in line:'//TRIM(ADJUSTL(LineStr))//'.'
      call MESSAGES(214,ErrMsg1,'fatal',INFO(1),INFO(2))   
	  return

      !test for existence of data file
201   rewind(u(m)%lu)
      read(u(m)%lu,'(A1)', err=770)dirbuch
      rewind(u(m)%lu)

      select case (u(m)%imode)
        case (0:1)                  !USER DEFINED DATA / WEATHER DATA
          ifile(1) = 0
!A MAXIMUM  OF maxcom COMMENT LINES IS ALLOWED; COUNT NUMBER
          rewind(u(m)%lu)
          do i=1,maxcom
            read(u(m)%lu,*,end=200) cdum
            read(cdum,'(i1)',err=100,end=200) ix
            if (xflag==0) xflag=1
            cycle
  100       ifile(1)=ifile(1)+1
            if (xflag==1) goto 801     !DATA LINE FOLLOWED BY COMMENT LINE
          enddo
  200     if (ifile(1)>maxcom) goto 802  !MORE COMMENTS THAN ALLOWED
          rewind(u(m)%lu)


        case (2)                  !TMY2-DATA
         read(u(m)%lu,'(33X,i3,x,a,x,i2,1X,i2,1X,A1,1X,i3,1X,i2)',err =772) ish,dirLat,ilatd,ilatm,dirbuch,ilongd,ilongm
          r(m)%latit   = ilatd+ilatm/60.
          if ( (dirLat == 's') .or. (dirLat == 'S') ) then
            r(m)%latit = -r(m)%latit
          endif
          if (dirbuch=="W") then
            r(m)%shift = -(ilongd+ilongm/60.)-ish*15.
          else
            r(m)%shift =  (ilongd+ilongm/60.)-ish*15.
          endif
          ifile(1) = 1

       case (3)                  !TRY-DATA
         do i=1,2; read(u(m)%lu,*); enddo
         read(u(m)%lu,'(A5,f6.2,13X,f6.2)',err=773) dirword, r(m)%latit, r(m)%shift
         if (dirword.ne.'LAGE:') goto 772
         r(m)%shift = -15.+r(m)%shift
         ifile(1) = 24

       case (9:)              !VDI-DATA
         read(u(m)%lu,'(29X,A5)',err=779) dirword
         if (dirword.ne.'T=3.7') goto 779
         rewind(u(m)%lu)
         ifile(1) = (designmonth-1)*26+2
         ifile(2) = 0
         ifile(3) = 24
         ndata    = 24
         goto 700

       case default
      end select

      rewind(u(m)%lu)

      do i=1,ifile(1); read(u(m)%lu,*); enddo

!NOW TRY TO SKIP THE DATALINES FROM T1 (TIME OF FIRST DATA LINE) TO
!NSKIP (DATA LINE AT TIME time0, which is 1 o'clock of first simulation day) AND COUNT SUCCESSFULL TRIES
      if (nskip>0) then !SIMULATION BEGINS LATER THAN DATA FILE
       do i=1,nskip
        read(u(m)%lu,'(A)',end=300,err=780) cdum
        if (len_trim(cdum)==0) then
          WRITE(LineStr,*) i
          ErrMsg1 = 'Error ocurred in line:'//TRIM(ADJUSTL(LineStr))//'.'
          call MESSAGES(214,ErrMsg1,'fatal',INFO(1),INFO(2))   
		  return
		endif
        read(cdum,'(e2.0)',end=300,err=780) xtest
	    ifile(2)=ifile(2)+1
       enddo
      endif

!NOW TRY TO READ IN ndata DATA LINES AND COUNT SUCCESSFULL TRIES
  300 do i=1, ndata
        read(u(m)%lu,'(A)',end=400,err=780) cdum
        if (len_trim(cdum)==0) then
          WRITE(LineStr,*) i
          ErrMsg1 = 'Error ocurred in line:'//TRIM(ADJUSTL(LineStr))//'.'
          call MESSAGES(214,ErrMsg1,'fatal',INFO(1),INFO(2))   
		  return
        endif

        read(cdum,'(e2.0)',end=400,err=780) xtest
        ifile(3)=ifile(3)+1
      enddo

!FOR SPLINE INTERPOLATION TWO MORE LINES ARE NEEDED
  400 do i=1,2
        read(u(m)%lu,'(f2.0)',end=500,err=780)xtest;ifile(5) = ifile(5)+1
      enddo

  500 if (u(m)%ixh==1 .and. ndata/=ifile(3)) then        !COUNT ALL DATA LINES
      do
        read(u(m)%lu,*,end=600); ifile(5) = ifile(5)+1
      enddo
      endif
  600 if (u(m)%ixh==0 .and. ifile(3)<ndata) goto 803			!no repeat allowed but not enough data

      ndata=ifile(2)+ifile(3)									!NUMBER OF DATA LINES
      ttot= u(m)%ddelt*ndata									!DURATION OF DATA ARRAY

      if (readerr<1 .and. u(m)%irad>0 .and. ifile(3) <ifile(4).and.u(m)%ixh==1) then  !DOES REPEAT READ FIT?
        if (mod(ttot,24.)>0.) goto 805                 !DATA ARE NOT MULTIPLE OF 24 H
        ifile(2)=0                                      ! NO SKIP
        ifile(3)=ndata                                  ! NO SKIP
        ifile(4)=ndata                                  ! NO SKIP  !TL TW commented out 24/2/05 bug repeat 2 days dont know if correct?
        ifile(5)=0                                      ! NO SKIP
      endif

  700 rewind(u(m)%lu)
      if (ifile(3)<ifile(4).and.ifile(2)>0) then        !SIMULATION STARTS LATER THAN DATA AND NEED REWIND
        nskip = 0; xflag=0
      endif
      r(m)%latit = r(m)%latit*rdconv
      r(m)%latit = sign(max(abs(r(m)%latit),1.e-6),r(m)%latit)
      u(m)%ifile = max(u(m)%ifile,ifile)

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('infiletest')
      if (time_report) call write_clock_time ('infiletest')
! ---------------------------------------------------------------------

      end subroutine infiletest
!*****************************************************************************


!*****************************************************************************
      subroutine dataread
!*****************************************************************************

! READ THE DATA FROM FILE TO ARRAY.

      integer ::  cc
	  character :: LineStr
	  character(len=200) :: ErrMsg1

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('dataread')
      endif
! ---------------------------------------------------------------------

! THIS PART IS LENGTHY; BECAUSE THE SPLINE INTERPOLATION NEEDS TWO MORE DATA
      cc =  0                                            !AT BEGIN AND END OF DATA ARRAY - SORRY
      do i=1,ifile(1);read(u(m)%lu,*);cc=cc+1;enddo      !SKIP COMMENT LINES
      if (u(m)%ixh==0) then                              !NO REWIND  u(m)%ixh
        do i=1,ifile(2)-2;read(u(m)%lu,*);cc=cc+1;enddo  !SKIP DATA LINES

        cc=cc+1
        if (u(m)%iform==0) then
          read(u(m)%lu,*,err=800)(u(m)%p(1,u(m)%ro(j)),j=1,u(m)%nn)
        else
          read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(1,u(m)%ro(j)),j=1,u(m)%nn)
        endif
        if (ifile(2)<2)then; backspace(u(m)%lu);cc=cc-1;endif

        cc=cc+1
        if (u(m)%iform==0) then
          read(u(m)%lu,*,err=800)(u(m)%p(2,u(m)%ro(j)),j=1,u(m)%nn)
        else
          read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(2,u(m)%ro(j)),j=1,u(m)%nn)
        endif
        if (ifile(2)==0)then; backspace(u(m)%lu); cc=cc-1; endif

        do i=3,updim-2
          cc=cc+1
          if (u(m)%iform==0) then
            read(u(m)%lu,*,err=800)(u(m)%p(i,u(m)%ro(j)),j=1,u(m)%nn)
          else
            read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(i,u(m)%ro(j)),j=1,u(m)%nn)
          endif
        enddo
        if (ifile(5)<2)then; backspace(u(m)%lu);cc=cc-1;endif

        cc=cc+1
        if (u(m)%iform==0) then
          read(u(m)%lu, *,err=800) (u(m)%p(updim-1,u(m)%ro(j)),j=1,u(m)%nn)
        else
          read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(updim-1,u(m)%ro(j)),j=1,u(m)%nn)
        endif
        if (ifile(5)<1)then; backspace(u(m)%lu);cc=cc-1;endif

        cc=cc+1
        if (u(m)%iform==0) then
          read(u(m)%lu,*,err=800)(u(m)%p(updim,u(m)%ro(j)),j=1,u(m)%nn)
        else
          read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(updim,u(m)%ro(j)),j=1,u(m)%nn)
        endif
      else                                    ! REWIND OPTION
        ndata2=mod(nint((time0-u(m)%t1)/u(m)%ddelt),ndata)
        ndata2=max(ndata2,0)  !Problem with START=0
        ndata1=ndata-ndata2
        ifile(2)=max(1,ceiling((time0-u(m)%t1)/u(m)%ddelt/ndata))*ndata-ndata1
        do i=3+ndata1,3+ndata-1
          if (u(m)%iform==0)then
          read(u(m)%lu,  *,err=800)(u(m)%p(i,u(m)%ro(j)),j=1,u(m)%nn)
          else
            read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(i,u(m)%ro(j)),j=1,u(m)%nn)
        endif
        enddo

        do i=3,3+ndata1-1
          if (u(m)%iform==0) then
          read(u(m)%lu,*,err=800)(u(m)%p(i,u(m)%ro(j)),j=1,u(m)%nn)
		else
            read(u(m)%lu,u(m)%stm,err=800) (u(m)%p(i,u(m)%ro(j)),j=1,u(m)%nn)
          endif
        enddo
        u(m)%p(1,:)=u(m)%p(updim-3,:); u(m)%p(2,:)=u(m)%p(updim-2,:)
        u(m)%p(updim-1,:)=u(m)%p(3,:); u(m)%p(updim,:)=u(m)%p(4,:)
      endif

! DO UNIT CONVERSIONS FOR INITIAL DATA AND SET TIME VALUES
      u(m)%t1 = u(m)%t1+ifile(2)*u(m)%ddelt
      do i=1,updim
        u(m)%p(i,-1) = u(m)%t1+(i-3.+.5*u(m)%tshift)*u(m)%ddelt
        do j=1,u(m)%nn
          u(m)%p(i,j)=u(m)%prodf(j)*u(m)%p(i,j)+u(m)%addf(j)
        enddo
      enddo
      u(m)%p(2,-1)       = u(m)%p(3,-1)-u(m)%ddelt
      u(m)%p(1,-1)       = u(m)%p(2,-1)-u(m)%ddelt
      u(m)%p(updim-1,-1) = u(m)%p(updim-2,-1)+u(m)%ddelt
      u(m)%p(updim  ,-1) = u(m)%p(updim-1,-1)+u(m)%ddelt
      u(m)%p(:,0)=u(m)%p(:,-1)

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('dataread')
      if (time_report) call write_clock_time ('dataread')
! ---------------------------------------------------------------------

      return

800   WRITE(LineStr,*) i
      ErrMsg1 = 'Error ocurred in line:'//TRIM(ADJUSTL(LineStr))//'.'
      call MESSAGES(214,ErrMsg1,'fatal',INFO(1),INFO(2))   
	  return


      end subroutine dataread
!*****************************************************************************


!*****************************************************************************
      subroutine update_time_sun_radiation
!*****************************************************************************

      real :: t2a,t2,t2e,t3,ddelt,tia,tie,ydum
	integer irow
	  real :: xlastinterpolation(MAXU)   ! defined xlastinterpolation(m) : 25.10.2005  YY

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('update_time_sun_radiation')
      endif
! ---------------------------------------------------------------------

! UPDATE TIME FOR REPEATED DATA
       updim = u(m)%nupdim

      if (xtime>u(m)%p(updim-2,-1)+0.5*u(m)%ddelt .and. xtime<tfinal+0.5*u(m)%tshift*u(m)%ddelt .and. u(m)%ixh>0) then
        timeadd = u(m)%p(updim-2,-1)-u(m)%p(2,-1)
        u(m)%p(:,-1)=u(m)%p(:,-1)+timeadd
        u(m)%p(:,0) =u(m)%p(:,-1)
        u(m)%sunsetflag=0; n_low(m,:)=2
      endif

! TIME VALUES FOR RADIATION HAVE TO ACCOUNT FOR SUNRISE AND SUNSET
! ONLY WHEN TIME ARRAY CHANGES AND RADIATION INTERPOLATION CHOSEN
      if (u(m)%sunsetflag==0.and.u(m)%irad>=1) then
        u(m)%sunsetflag=1; timep=3
! TL 30042002        
! time for definition of day changed from p(3,-1) to p(4,-1)
! before: day= mod(ceiling((u(m)%p(3,-1)+0.001)/24.),365)
        day= mod(ceiling((u(m)%p(4,-1)+0.001)/24.),365)
        if (day<1) day=365
        do while(timep<=updim-2)
          declin = 0.40928*sin((284.+day-1)*piconv)
          if (u(m)%imode>9) declin = .40928*sin((284.+r(m)%designday)*piconv)       !VDI-Mode
          tanprod   = max(-1.,min(1.,-tan(declin)*tan(r(m)%latit)))
          riseangle = acos(tanprod)/rdconv
          sunrise =-(riseangle+r(m)%shift)/15.-(1-u(m)%ias)*ecc(mod(u(m)%p(timep,-1),8760.))+24.*day -12.+floor(xtime/8760.)*8760
          sunset  = sunrise+riseangle/7.5
          do jj=1,3; if (u(m)%p(jj,0)>sunrise) u(m)%p(jj,-2)=1; enddo
          ilocb=timep

          if (u(m)%p(timep,0)<sunrise) then
            do while ((u(m)%p(timep,0)-sunrise)/u(m)%ddelt<-1. .and. timep<=updim-1)
              u(m)%p(timep,-2)=0; timep=timep+1
            enddo

            if (abs(u(m)%p(timep,0)-sunrise)<u(m)%ddelt) then
              u(m)%p(timep,0) = sunrise
              u(m)%p(timep,-2) = 2
            endif
            if (timep<=updim-2) timep=timep+1
          end if

          if (u(m)%p(timep,0)<sunset) then
            do while((u(m)%p(timep,0)-sunset)/u(m)%ddelt<1. .and. timep<=updim-2)
              u(m)%p(timep,-2)=1; timep=timep+1
            enddo
            u(m)%p(timep-1,0) = sunset
            u(m)%p(timep-1,-2) = 3
          end if

          iday= mod(u(m)%p(timep,-1)/24.,365.)
          do while(mod(u(m)%p(timep,-1)/24.,365.)==iday .and. timep<updim-1)
            timep=timep+1
          enddo
          iloce=timep-1
          do i=ilocb,iloce
            do j=1,u(m)%nn           !CHECK RADIATION DATA FOR VALUES AT NIGHT
              if (u(m)%p(i,-2)/=1..and.u(m)%ipol(j)>2.and.u(m)%p(i,j)>3.6 .and. u(m)%p(i,0)>u(m)%p(i,-1)+.5*u(m)%ddelt .and. u(m)%p(i,0)<tfinal) then
                if (falseradflag(m)<11) iwarn = iwarn + 1
                falseradflag(m) = falseradflag(m)+1
                shloc = r(m)%shift/15.+(1-r(m)%ias)*ecc(u(m)%p(i,-1))
                if (falseradflag(m)<11) then
!CMH_09_05_06                  write(luw,1550) falseradflag(m), info(1),info(2),u(m)%p(i,-1)+.5*u(m)%ddelt,u(m)%p(i,j), j, sunrise,sunset,shloc
                  write(luw,1550) falseradflag(m), info(1),info(2),u(m)%p(i,-1)+.5*u(m)%ddelt,u(m)%p(i,j), sunrise,sunset,shloc
	          endif
              endif
           enddo
          enddo
          timep=timep+1; day = day+1
        enddo
        if (falseradflag(m)>0) write(luw,1560) falseradflag(m)
      endif

! THE FOLLOWING PART IS EXECUTED AT THE BEGINNING OF EVERY NEW DAY; FOR EVERY
! TIME STEP OF THE DAY RADIATION IS INTERPOLATED AND SCALED, THEN DATA HAVE
! JUST TO BE READ FROM THE ARRAY DURING SIMULATION TIME OF THIS DAY
! changed xlastinterpolation to xlastinterpolation(m) : 25.10.2005  YY
      if (info(7)==-1 .or.( nint(mod(xtime,24.)/delt)==1 .and. xlastinterpolation(m) /= xtime)) then !FIRST CALL OR NEW DAY
		xlastinterpolation(m) = xtime              
        iradcol=0; n_lo = n_low(m,-1)
        do i=1,u(m)%nn; if (u(m)%ipol(i)==3) then  !INTEGRATE FOR RADIATION INTERPOLATION
          lastdaycut=.false.
          frstdaycut=.false.
          nfitloc=u(m)%nfit-nint(mod(xtime,24.)/u(m)%ddelt) !SUBSTRACT BEGINTIME LATER THAN BEGIN OF DAY
          nlastday=nint((tfinal-xtime)/u(m)%ddelt)+1
          if (nlastday<nfitloc) then             !SIMUL STOPS BEFORE END OF DAY
            nfitloc=nlastday; lastdaycut=.true.
          endif
          iradcol=iradcol+1
          if (iradcol>u(m)%irad) return
          do k=1,nday; u(m)%dp(k,iradcol)=0.; enddo
          do while (u(m)%p(n_lo+1,0)<=xtime); n_lo=n_lo+1; enddo
          if (u(m)%p(n_lo,-2)>0.) frstdaycut=.true.

          u(m)%xfit=0.; u(m)%yfit=0.              !Initialize

          do k = 1,min(nfitloc,n_lo+nfitloc-1)
	      if (u(m)%p(n_lo+k-1,i)<0.) then !test for negative radiation values
              ydum = u(m)%p(n_lo+k-1,i)
              u(m)%p(n_lo+k-1,i) = 0.       !if negative rad-value=0
              irow = n_lo+k-2+u(m)%ifile(2)
			  call MESSAGES(-1,'Radiation data read from user-defined weather data were set to zero due to zero value.','warning',INFO(1),INFO(2))
            endif
            if (u(m)%p(n_lo+k-1,i)>0..and.u(m)%p(n_lo+k-1,-2)<1.) then !test for  rad at night
              ydum = u(m)%p(n_lo+k-1,i)
              u(m)%p(n_lo+k-1,i) = 0.       !if no sun rad-value=0
              irow = n_lo+k-2+u(m)%ifile(2)
              write(luw,1571) info(1),info(2),trim(u(m)%ofile),irow,ydum
            endif
            u(m)%xfit(k)=u(m)%p(n_lo+k-1,-1)
            u(m)%yfit(k)=u(m)%p(n_lo+k-1,i)
            if (u(m)%p(n_lo+k-1,-2)>2.) then
              suns = u(m)%p(n_lo+k-1,0)
              sunsmod=suns-mod(suns,delt)+delt
            else if (u(m)%p(n_lo+k-1,-2)>1.) then
              sunr = u(m)%p(n_lo+k-1,0)
              sunrmod=sunr-mod(sunr,delt)
            endif

          enddo
          if (suns<xtime.and.lastdaycut) suns=tfinal

!WORK TROUGH ALL TIME INTERVALS OF THIS DAY, FROM RAD=0 TO RAD=0
          nradi=0
          do while (nfitloc>0)
            if (nradi>0) nradi=nradi-1
            u(m)%xfit=cshift(u(m)%xfit,nradi)
            u(m)%yfit=cshift(u(m)%yfit,nradi)
            nfitloc=nfitloc-nradi
            if (.not. u(m)%yfit(1)>0) then     !IF SIMULATION STARTS WHILE RAD ON
              do while (.not.(u(m)%yfit(1)==0. .and. u(m)%yfit(2)/=0.).and.nfitloc>0) ! THIS PART WOULD NOT WORK
                u(m)%xfit=cshift(u(m)%xfit,1)
                u(m)%yfit=cshift(u(m)%yfit,1)   !SEEK RAD. != 0
                nfitloc=nfitloc-1
              enddo
            endif


            if (.not.nfitloc>0) exit
            nradi = 1
            do while (u(m)%yfit(nradi+1)>0..and.nfitloc>nradi)
              nradi=nradi+1        !SEEK INTERVAL WHERE RADIATION VALUES > 0
            enddo
            if (u(m)%xfit(nradi+1)<=tfinal.and.nradi<nfitloc) nradi = nradi+1
            if (lastdaycut) nfitloc=0

            kbeg = max(1,ceiling(mod(u(m)%xfit(1)  ,24.)/delt))
            kend = max(1,ceiling(mod(u(m)%xfit(nradi),24.)/delt))
            timeloc = xtime-nint(mod(xtime,24.)/delt)*delt
            if (kbeg>=kend) exit

            realsum=0.                  !CALCULATE SUM OF GIVEN RADIATION DATA
            do kkh=1,nradi
             realsum=realsum+u(m)%yfit(kkh)*u(m)%ddelt
            enddo
            if (frstdaycut) then
              realsum=realsum-u(m)%yfit(1)*u(m)%ddelt/2.
              frstdaycut=.false.
            endif
            if (lastdaycut.and.u(m)%xfit(nradi)>=tfinal) then
              realsum=realsum-u(m)%yfit(nradi)*u(m)%ddelt/2.
            endif

!correct for sunrise and sunset
            xfitloc=u(m)%xfit; yfitloc=u(m)%yfit
            ddelt = u(m)%ddelt
               t2 = u(m)%xfit(2)
              t2a = t2  - ddelt/2
              t2e = t2a + ddelt

            if (sunr>t2a .and. sunr<t2e) then
	          yfitloc(2) = yfitloc(2)*ddelt/(t2e-sunr)
              t2a = sunr
              t2  = (t2a+t2e)/2.
              xfitloc(2) = t2
              !data are not correct? if to high: restrict to 60% of following value
              if (yfitloc(2) > 0.6*yfitloc(3)) yfitloc(2) = 0.6*yfitloc(3)
            end if

            if (nradi>2) then
               tn1 = xfitloc(nradi-2)
			else
			  tni = xfitloc(1)
			endif

               tn2 = xfitloc(nradi-1)
              tn2a = tn1  + ddelt/2
              tn2e = tn2a + ddelt

            if (sunsmod>tn2a .and. sunsmod<tn2e) then
              yfitloc(nradi-1) = yfitloc(nradi-1)*ddelt/(sunsmod-tn2a)
              tn2e = sunsmod
              tn2  = (tn2a+tn2e)/2.
              xfitloc(nradi-1) = tn2
              !data are not correct?
              if (yfitloc(nradi-1) > 0.6*yfitloc(nradi-2)) yfitloc(nradi-1) = 0.6*yfitloc(nradi-2)
            end if

!only one or two data points occur?
            xfitloc(1) = t2a
            xfitloc(nradi)= tn2e


            do idp = kbeg, kend
              xtemp = timeloc+delt*(idp-.5)
              tia = xtemp - delt/2.
              tie = xtemp + delt/2.

              if (tie <=  t2a) cycle
              if (tia >= tn2e) cycle

              if (tia < t2a .and. tie > t2a) then     !part of rad in sim interval at rad start
                  call akimaip(luw,nradi,xfitloc,yfitloc,1,0.5*(t2a+tie),radv(1))
                  ydum = max(radv(1),0.)
                  u(m)%dp(idp,iradcol) = ydum*(tie-t2a) / delt
              end if

              if (tia < tn2e .and. tie > tn2e) then !part of rad in sim interval at rad end
                  call akimaip(luw,nradi,xfitloc,yfitloc,1,0.5*(tia+tn2e),radv(1))
                  ydum = max(radv(1),0.)
                  u(m)%dp(idp,iradcol) = ydum*(tn2e-tia) / delt
              endif

              if (tia >= t2a .and. tie <= tn2e) then !everything quite normal case
                call akimaip(luw,nradi,xfitloc,yfitloc,1,xtemp,radv(1))
                u(m)%dp(idp,iradcol)=max(radv(1),0.)
              endif
            enddo

            tempsum=0.
!tl, may 2002: counter altered for small radiation error when radiation goes down to 0 during daytime            
!           do ic = kbeg,kend
            do ic = kbeg+1,kend
              tempsum=tempsum+u(m)%dp(ic,iradcol)*delt
            enddo
            if (tempsum==0.) tempsum=1.; factkorr=realsum/tempsum
            do ic = kbeg,kend
              u(m)%dp(ic,iradcol)=u(m)%dp(ic,iradcol)*factkorr
            enddo
          enddo
        endif; enddo
      endif


! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('update_time_sun_radiation')
      if (time_report) call write_clock_time ('update_time_sun_radiation')
! ---------------------------------------------------------------------

      return



1550  FORMAT(//' ***** WARNING ',i2,' *****  FROM UNIT ',I3,' TYPE ',I3,' DATA READER',/,' RADIATION AT TIME',F8.2,' HAS A VALUE OF ',F8.2,' BUT SUN IS DOWN',/,' -- SUNRISE AT ',F8.2,', SUNSET AT ',F8.2,/,' (SOLAR TIME = TIME + ',F8.2,')')
1560  FORMAT(//' Radiation while sun is down found ',i4,' times',/,' only the first 10 warnings with data info are printed ',/,' check: - Time of first data line',/,'        - Longitude and Latitude',/,' REMARK: these radiation values are set to ZERO !')
1570  FORMAT(//' ***** WARNING FROM UNIT ',I3,' TYPE ',I3,' DATA READER',/,' RADIATION DATA READ IN FROM FILE ==>',A,/,' AT LINE #',I5,' ARE SET TO ZERO DUE TO NEGATIVE VALUE',F8.2,' W/m²')
1571  FORMAT(//' ***** WARNING FROM UNIT ',I3,' TYPE ',I3,' DATA READER',/,' RADIATION DATA READ IN FROM FILE ==>',A,/,' AT LINE #',I5,' ARE SET TO ZERO SINCE SUN IS DOWN, ORIGINAL ','VALUE:',F8.2,' W/m²')

      end subroutine update_time_sun_radiation
!*****************************************************************************

!*****************************************************************************
      END SUBROUTINE TYPE109
!*****************************************************************************


!*****************************************************************************
      subroutine akimaip(iu,l0,x,y,n0,u,v)
!*****************************************************************************
! Splineinterpolation nach Akima
! copied and recycled by Thomas Lechner, Transsolar, 15-2-96
! changed to radiation interpolation: fitted values at beginning and end of
! data interval don't swing below 0
! Source:  Interpolation and Smooth Curve Fitting Based on Local Procedures
!          Hiroshi Akima
!          Communications of the ACM
!          October 1972, Vol. 15, No. 10, pp. 914-918
! Input data:
!       iu      standard output unit
!       l       number of input values, >=2
!       x       Array of length l, x-Werte input
!       y       Array of length l, y-Werte input
!       n       number of x-values interpolation will be done for
!       u       Array of length n, x-values output
! Output data:
!       v       Array of length n, y-values output

! Export this subroutine for its use in external DLLs
!dec$attributes dllexport :: akimaip

! declarations
! type109_time_report: 24.02.2005 TW ------------------------------
      use type109_time_report

integer :: I, J, K, L0, N0, IPV,IU
  
real(4) :: Y5, A4, A1, A5, W2, W3, SW, T4, SA, T3, W4 
real(4) :: UK, Q2, Q3, DX, Q0, Q1, X3, Y3, X4, Y4, A3, X2, Y2, A2, X5 
real(4) :: m1,m2,m3,m4,m5

real(4) :: x(l0),y(l0),u(n0),v(n0)

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('akimaip')
      endif
! ---------------------------------------------------------------------

! initialise:
      if (l0<2) goto 2000
      if (n0.le.0) goto 2010
      do i=2,l0; if (x(i-1).ge.x(i)) goto 2020; enddo
      ipv = 0

! main loop
      do k=1,n0
      uk = u(k)
      if      (l0 == 2    ) then; i = 2
      else if (uk >= x(l0)) then; i = l0+1
      else if (uk <  x(1) ) then; i = 1
      else
        i = l0; do while (uk<x(i-1)); i=i-1; enddo
      endif

! check if i = ipv
      if (i.eq.ipv) goto 70
      ipv = i

! pick up x,y and estimate if necessary
      j = i; if (j.eq.1) j=2; if (j.eq.l0+1) j=l0
      x3 = x(j-1); y3 = y(j-1); x4 = x(j); y4 = y(j)
      a3 = x4-x3; m3 = (y4-y3)/a3

      if (l0==2) then                 !linear interpolation
        m2 = m3; m4 = m3; m1 = m2+m2-m3; goto 50
      endif
      if (j==l0) then                 !last x-value, m4 set to 0
        x2 = x(j-2); y2 = y(j-2); a2 = x3-x2
        m2 = (y3-y2)/a2; m4=0; goto 50
      endif
      if (j==2) then                  !first value, m2 set to 0
        x5 = x(j+1); y5 = y(j+1); a4 = x5-x4
        m4 = (y5-y4)/a4; m2=0; goto 40
      endif

      x2 = x(j-2); y2 = y(j-2); a2 = x3-x2
      x5 = x(j+1); y5 = y(j+1); a4 = x5-x4
      m2 = (y3-y2)/a2; m4 = (y5-y4)/a4

   40   if (j>3) then
        a1 = x2-x(j-3); m1 = (y2-y(j-3))/a1
      else
        m1 = m2+m2-m3
      endif
      if (j.ge.l0-1) then
        m5 = m4+m4-m3
      else
        a5 = x(j+2)-x5; m5=(y(j+2)-y5)/a5
      endif

! numerical differentiation
   50   if (i==1) then
        w2 =abs(m4-m3); w3 = abs(m2-m1); sw = w2+w3
        if (sw==0.0) then; w2 = .5; w3 = .5; sw = 1.; endif
        t4 = (w2*m2+w3*m3)/sw; sa = a3+a4
        t3 = .5*(m1+m2-a4*(a3-a4)*(m3-m4)/sa/sa)
        x3 = x3-a4; y3 = y3-m2*a4; a3 = a4; m3 = m2
      else if (i==l0+1) then
        w3 = abs(m5-m4); w4 = abs(m3-m2); sw = w3+w4
        if (sw==0.0) then; w3 = .5; w4 = .5; sw = 1.; endif
        t3 = (w3*m3+w4*m4)/sw; sa = a2+a3
        t4 = .5*(m4+m5-a2*(a2-a3)*(m2-m3)/sa/sa)
        x3 = x4; y3 = y4; a3 = a2; m3 = m4
      else
        w2 =abs(m4-m3); w3 = abs(m2-m1); sw = w2+w3
        if (sw==0.0) then; w2 = .5; w3 = .5; sw = 1.; endif
        t3 = (w2*m2+w3*m3)/sw
        w3 = abs(m5-m4); w4 = abs(m3-m2); sw = w3+w4
        if (sw==0.0) then; w3 = .5; w4 = .5; sw = 1.; endif
        t4 = (w3*m3+w4*m4)/sw
      endif

!determination of coefficients
      q2 = (2.*(m3-t3)+m3-t4)/a3
      q3 = (-m3-m3+t3+t4)/a3**2

! computation of the polinominal
   70 dx = uk-x3
      q0 = y3; q1 = t3
      v(k) = q0+dx*(q1+dx*(q2+dx*q3))
      enddo

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('akimaip')
      if (time_report) call write_clock_time ('akimaip')
! ---------------------------------------------------------------------

      return
 2000 write(iu,*)'weniger als 2 Eingabepunkte in sub AKIMAIP'; return
 2010 write(iu,*)'weniger als 1 Ausgabepunkt in sub AKIMAIP'; return
 2020 write(iu,*)'***** possibly interpolation problem at time',x(i),'***** (internal information)'; return
      end subroutine akimaip




      SUBROUTINE rad_proc(time,rin,rres)
! special changes of original radiation processor:
! using fortran 90 syntax, special repeat of a design day for VDI 2078
! Thomas Lechner, Transsolar February 2001
!
! 
!  TL, 27.05.02   peak at the end of a day -fixed -test ok
!   
!*******************************************************************************
!   
!   THIS ROUTINE FINDS THE POSITION OF THE SUN IN THE SKY AND
!   ESTIMATES THE RADIATION INCIDENT ON AN ARBITRARY NUMBER OF SURFACES OF
!   ANY ORIENTATION.
!
!   THIS ROUTINE WAS ORIGINALLY WRITTEN IN THE DARK AGES OF TRNSYS HISTORY.  
!   IT WAS EXPANDED IN EARLY 1989 BY DOUG REINDL AND 'CLEANED' BY RUTH URBAN.
!   2/92 ERROR FIX -- JWT
!   2/93 RADIATION SMOOTHING -- RAG
!   2/93 PREPARED FOR INCLUSION IN 14.0 -- JWT
!  
!   THIS ROUTINE IS CALLED BY: TYPE 169-DATAREADER BY THOMAS LECHNER
!
!*******************************************************************************
use TrnsysFunctions ! Use TRNSYS global functions 
use TrnsysData, ONLY : LUR,LUW,IFORM,LUK ! Use TRNSYS global constants 

      use Type109RadiationData
! type109_time_report: 24.02.2005 TW ------------------------------
      use type109_time_report
      type (rad_proc_call), intent(inout) :: rin
      type (rad_proc_call)                :: r

      integer :: day          !CURRENT DAY
      integer :: ip           !POINTER TO FIRST INPUT FOR A GIVEN SURFACE
      integer :: iwarn        !ERROR COUNTER
      integer :: jp           !POINTER TO FIRST OUTPUT FOR A GIVEN SURFACE
      integer :: nbin         !INTEGER VARIABLE
      integer :: np           !NUMBER OF PARAMETERS          
      integer :: norad        !NO RADIATION-FLAG
      integer :: i,is         !LOOP COUNTER

      real :: a1,ai       !VARIABLE
      real :: alf         !SOLAR AZIMUTH - SURFACE AZIMUTH
      real :: angle       !FUNCTION TO CALCULATE SOLAR TIME HOUR ANGLE SEE BELOW
      real :: aux         !VARIABLE
      real :: axazm       !AZIMUTH OF CURRENT SURFACE OR TRACKING AXIS
      real :: axslp       !SLOPE OF CURRENT SURFACE OR TRACKING AXIS
      real :: azm         !ACTUAL AZIMUTH OF CURRENT SURFACE
      real :: b1          !VARIABLE
      real :: cc          !PRODUCT OF COS(LATITUDE)*COS(DECLINATION)
      real :: cosasl,cosdec,coshr,cosslp,costtp,costt,coszen !VARIABLES
      real :: cwew        !VARIABLE
      real :: decl        !DECLINATION
      real :: dgconv=57.2958  !CONVERT RADIANS/DEGREE
      real :: delt        ! TIME STEP DURATION (HOURS)
      real :: dfract      !DIFFUSE FRACTION (HD/HHOR)
      real :: ecc         !ACCOUNTS FOR THE VARIATION OF EXTRATERRSTRIAL RADIA-
                          !TION DUE TO VARIATION IN SUN TO EARTH DISTANCE
      real :: eps,epsiln  !VARIABLE
      real :: et          !TIME SHIFT DUE TO ECCENTRICITY OF THE EARTH'S ORBIT
      real :: eot         !FUNCTION TO CALCULATE SOLAR TIME HOUR ANGLE SEE BELOW
      real :: f           !VARIABLE
      real :: gam         !AZM - AXAZM
      real :: hb          !BEAM HORIZONTAL RADIATION
      real :: hbeam       !BEAM RADIATION ON TILTED SURFACE
      real :: hd          !DIFFUSE HORIZONTAL RADIATION
      real :: hdiff       !DIFFUSE RADIATION ON TILTED SURFACE
      real :: hdn         !DIRECT NORMAL RADIATION
      real :: hextra,hx   !HORIZONTAL EXTRA-TERRESTRIAL RADIATION
      real :: hgrf        !ISOTROPIC GROUND REFL. RADIATION ON A TILTED SURFACE
      real :: hhor        !TOTAL HORIZONTAL RADIATION
      real, dimension(8)      :: p11,p12,p13,p21,p22,p23 !USED IN PEREZ SKY MODEL
      real :: p1,p2       !VARIABLE
      real :: pi=3.1415927
      real :: piconv=.0172028 !2*PI/365.242
      real :: rdconv=.0174533 !CONVERT DEGREE/RADIANS
      real :: rb          !VARIABLE
      real :: rh          !RELATIVE HUMIDITY
      real :: sazm        !SOLAR AZIMUTH ANGLE
      real :: sc=4871.15  !SOLAR CONSTANT
      real :: scube       !VARIABLE
      real :: sinazm,sinasl,sindec,sinhr,sinzen,sinslp        !VARIABLES
      real :: skyb        !VARIABLE
      real :: slope       !ACTUAL SLOPE OF CURRENT SURFACE
      real :: ss          !SIN(LATITUDE) * SIN(DECLINATION)
      real :: ta          !VARIABLE
      real :: tanasl,tandec,tangam,tnprod,tanslp,tanzen   !VARIABLES
      real :: theta       !ANGLE OF INCIDENCE FOR BEAM RADIATION
      real :: time        !CURRENT SIMULATION TIME
      real :: tims,time0,tfinal  !FIXED TIME DUE TO DESIGN DAY (VDI)
      real :: wew         !VARIABLE
      real :: ws,w1,w2,w  !SUNSET/RISE HOUR ANGLE
      real :: x           !VARIABLE
      real :: xkt         !CLEARNESS INDEX
      real :: zenith      !ZENITH ANGLE OF SUN

      DOUBLE PRECISION, INTENT (OUT) :: rres(*)
      INTEGER*4 INFO(15)
      
!      COMMON /SIM/ TIME0,TFINAL,DELT,IWARN
!      COMMON /LUNITS/ LUR,LUW,IFORM,LUK

! DATA FOR PEREZ MODEL REPORTED IN SANDIA REPORT, 1988.
      DATA P11 /-0.196,0.236,0.454,0.866,1.026,0.978,0.748,0.318/
      DATA P12 /1.084,0.519,0.321,-0.381,-0.711,-0.986,-0.913,-0.757/
      DATA P13 /-0.006,-0.18,-0.255,-0.375,-0.426,-0.35,-0.236,0.103/
      DATA P21 /-0.114,-0.011,0.072,0.203,0.273,0.28,0.173,0.062/
      DATA P22 /0.18,0.02,-0.098,-0.403,-0.602,-0.915,-1.045,-1.698/
      DATA P23 /-0.019,-0.038,-0.046,-0.049,-0.061,-0.024,0.065,0.236/


!STATEMENT FUNCTIONS TO CALCULATE SOLAR TIME AND HOUR ANGLES. EQUATION OF TIME 
! FROM D.L. SIEBERS, TECH. REPORT  ME-HTL-75-2, PURDUE UNIVERSITY, LAFAYETTE
      eot(x) = -.1236*sin(x)+.004289*cos(x)-.1539*sin(2.*x)-.06078*cos(2.*x)
      angle(x) = ((mod(x+et,24.)-12.)*15.+r%shift)*rdconv

        delt  = sngl(getSimulationTimeStep())
        time0 = sngl(getSimulationStartTime())
        tfinal= sngl(getSimulationStopTime())

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('rad_proc')
      endif
! ---------------------------------------------------------------------

      r = rin
      if (r%ivdi==1) then
          tims = (r%designday-1)*24.+mod(time,24.)
        else
          tims = mod(time,8760.) !cmh_29_11_05
      endif



!CALCULATIONS TO BE DONE AT EVERY TIMESTEP
      call calc_every_step(norad)

      if (norad>0) then; call no_radiation; goto 300; endif   !RETURN
!      if (hextra<1.e-6.or.(r%radmode.ne.6.and.r%iges.le.0.) .or. (r%radmode==6.and.r%id.le.0.)) then      
      if ((hextra<1.e-6).or.(r%iges.lt.0.).or.(r%id.lt.0.)) then
        hextra = 0.; call no_radiation; goto 300              !NO RADIATION, RETURN
      end if

      call suncalc        !CALCULATE SUN POSITION      
      call beamdifftotal  !FIND RATE OF BEAM, DIFF AND TOTAL RAD ON HORIZONTAL

      !SET OUTPUTS FOR EXTRAT., SOLAR ZENITH, SOLAR AZIMUTH
      rres(1)=hextra; rres(2)=zenith*dgconv; rres(3) = sazm*dgconv
      rres(4) = hhor              !TOTAL RAD ON HORIZONTAL
      rres(5) = hb                !BEAM RAD ON HORIZONTAL
      rres(6) = hd                !DIFFUSE RAD ON HORIZONTAL
      rres(7) = 0.                !GROUND REFLECTED RAD ON HORIZONTAL
      rres(8) = rres(2)           !SOLAR ZNEITH = ANGLE OF INCIDENCE FOR HOR.
      rres(9) = 0.                !SOLAR AZIMUTH FOR HORIZONTAL

! CALCULATE BEAM RADIATION, TOTAL RADIATION AND INCIDENCE ANGLE FOR EACH SLOPE
      ip = 1; jp = 10

      loop_over_surfaces: do is = 1,r%nsurf
        call sloped_surfaces
          !TOTAL, BEAM AND DIFFUSE AND GROUND REFLECTED RADIATION, INCIDENCE ANGLE, SLOPE 
          !FOR FIRST SURFACE.
          rres(jp)  = hbeam+hdiff+hgrf
          rres(jp+1)= hbeam; rres(jp+2) = hdiff; rres(jp+3)=hgrf
          rres(jp+4)= theta; rres(jp+5) = slope
          jp = jp + 6
      enddo loop_over_surfaces


! PRINT COUNT OF TIMES BEAM RADIATION GREATER THAN EXTRATERRESTRIAL
300   if (time < tfinal-hdelt) return
      if (r%nerr/=0) then
          write(luw,312) r%infloc(1), r%infloc(2),r%nerr
      endif
      iwarn=iwarn+1; return

312   FORMAT(//2X,24H***** WARNING ***** UNIT,I3,6H TYPE ,I3/4X,47HBEAM RADIATION EXCEEDED EXTRATERRESTRIAL DURING,I7,10H TIMESTEPS)

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('rad_proc')
      if (time_report) call write_clock_time ('rad_proc')
! ---------------------------------------------------------------------
      
      contains

! *********************************************************************************
      subroutine beamdifftotal
! *********************************************************************************
! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('beamdifftotal')
      endif
! ---------------------------------------------------------------------

!FIND RATE OF BEAM,DIFF AND TOTAL RADIATION ON HORIZONTAL SURFACE.

      hhor = r%iges
      xkt = hhor/hextra

      radiation_mode: Select case (r%radmode)
        case (1) ! USE LIU AND JORDAN CORRELATION FOR BEAM AND DIFFUSE.
          xkt = min(xkt,.75)
          hd = hhor*(1.0045+((2.6313*xkt-3.5227)*xkt+.04349)*xkt)
          hd = max(hd,0.); hd = min(hhor,hd); hb = hhor - hd

        case (2) ! USE BOES CORRELATION FOR DIRECT NORMAL RADIATION
          hdn = max(0.,(1.3304*xkt-.3843)*sc)
          if (hdn>.739*sc) hdn = .739*sc
          hb = min(hdn*coszen,hhor); hd = hhor-hb

        case (3) ! USE ERBS CORRELATION FOR BEAM AND DIFFUSE
          xkt_mode3: select case (int(100.*xkt))
          case (:22);  hd = hhor*(1.0 - 0.09*xkt)
          case (80:);   hd = 0.165*hhor
          case default; hd = hhor*(.9511+xkt*(-.1604+xkt*(4.38+xkt*(-16.638+12.336*xkt))))
          end select xkt_mode3
          hb = hhor - hd

        case (4) ! USE REINDL REDUCED CORRELATION TO ESTIMATE HORIZONTAL DIFFUSE
          xkt_mode4: select case (int(100.*xkt))
          case (:30);  dfract = 1.020-.254*xkt+.0123*coszen
          case (78:);  dfract=.486*xkt-.182*coszen
          case default; dfract = 1.4-1.749*xkt+.177*coszen
          end select xkt_mode4 
          dfract = min(dfract,.97); dfract = max(dfract,.10)
          hd = dfract*hhor; hb = hhor-hd

        case (5) ! REINDL FULL CORRELATION (TA IN C, REL HUMIDITY AS FRACTION)
          ta = r%t_amb; rh = r%rh_amb*.01
          xkt_mode5: select case (int(100.*xkt))
          case (:30 ) 
            dfract = 1.-.232*xkt+.0239*coszen-.000682*ta+.0195*rh
          case (78:)
            dfract = .426*xkt-.256*coszen+.00349*ta+.0734*rh
          case default
            dfract = 1.329-1.716*xkt+0.267*coszen-.00357*ta+.106*rh
          end select xkt_mode5 
          dfract = min(dfract,.97); dfract = max(dfract,.1)
          hd = dfract*hhor; hb = hhor - hd
     
        case (6) ! INPUTS ARE BEAM AND DIFFUSE ON HORIZONTAL
          hb = r%ib; hd = r%id; hhor = hb + hd

        case (7) ! INPUTS ARE TOTAL (HORIZONTAL) AND DIRECT NORMAL
          hdn = r%ibdn; hb = max(hdn*coszen,0.)
          hd = max(hhor-hb,0.)

        case (8) ! INPUTS ARE TOTAL (HORIZONTAL) AND DIFFUSE (HORIZONTAL)
          hd = r%id; hb = hhor - hd

		case (9) ! INPUTS ARE BEAM NORMAL RADIATION (IBN) AND DIFFUSE (HORIZONTAL)
		  hdn  = r%ibdn		!Get data from file
		  hd   = r%id		!Get data from file
		  hb   = max(hdn*coszen,0.)
		  hhor = hb + hd

      end select radiation_mode

! do not allow negative radiation
      hb = max(hb,0.)
!cmh_05_12_05 bugfix 1129 - different radiation values in the 2nd year due to lack of accuracy
      if(hb.lt.0.00009) hb=0. 

! CHECK FOR BEAM RADIATION GREATER THAN EXTRATERRESTRIAL RADIATION
      if(hb > hextra .and. hb > 3.6) then
        if((tims-hdelt) .ge. time0) then
          r%nerr = r%nerr+1
          if (r%nerr.le.12.) then
            write (luw,202) r%infloc(1),r%infloc(2),hb,tims,hextra
            iwarn=iwarn+1
          endif
        endif
        hb = hextra
      endif
      
202   FORMAT(//2X,24H***** WARNING ***** UNIT,I3,5H TYPE,I3/4X,15HBEAM RADIATION:,F10.2,34H EXCEEDS EXTRATERRESTRIAL AT TIME ,F8.2,/4X,27HEXTRATERRESTRIAL RADIATION:,F10.2,5H USED)

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('beamdifftotal')
      if (time_report) call write_clock_time ('beamdifftotal')
! ---------------------------------------------------------------------

      end subroutine beamdifftotal
! ***************************************************************************


! *********************************************************************************
      subroutine calc_every_step(norad)
! *********************************************************************************
!CALCULATIONS TO BE DONE AT EVERY TIMESTEP
      integer norad   !NO RADIATION

 ! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('calc_every_step')
      endif
! ---------------------------------------------------------------------

      norad = 0
!DECLINATION AND OTHER DAY DEPENDENT QUANTITIES.
!cmh 29_11_05       day=ceiling(tims/24.); decl=.40928*sin((284.+day)*piconv) 
      day=MOD(ceiling(tims/24.),365); decl=.40928*sin((284.+day)*piconv) !insert MOD becuase decl differed for the 2nd year due to a lack of accuracy
      sindec=sin(decl); cosdec=cos(decl); tandec=tan(decl)
      tnprod=-tandec*r%tanlat; tnprod=min(tnprod,1.)
      tnprod=max(tnprod,-1.)
      ws=acos(tnprod); cc=r%coslat*cosdec; ss=r%sinlat*sindec
      ecc=(1.+(1-r%ias)*.033*cos(day*piconv))
      et=(1-r%ias)*eot((day-1.)*piconv)

!FIND HOUR ANGLES FOR START AND END OF TIMESTEP. THE PORTION
!OF THE TIMESTEP DURING WHICH THE SUN IS DOWN IS IGNORED.
      w1 = angle(tims-delt); w2 = angle(tims); cdelt=0    !data from reader correspond to time interval delt prceeding actual time
      if (.not.(w2>-ws.and.w1<ws)) then
        norad=1; goto 999 !NO SUN
      endif 
      w1 = max(w1,-ws); w2 = min(w2,ws)
!! tl 	if (abs(w2-w1)<cdelt) then  27 May 2002
	if (abs(w2-w1)<0.01*delt) then
        norad=1; goto 999 !NO SUN
      endif
      w  = (w1+w2)*0.5; coshr = cos(w); sinhr = sin(w)
      hx = cc*(sin(w2)-sin(w1))+ss*(w2-w1)
      hextra = sc*ecc*hx/(w2-w1)

! type109_time_report: 24.02.2005 TW ------------------------------
  999 if (time_report) call get_clock_on_exit('calc_every_step')
      if (time_report) call write_clock_time ('calc_every_step')
! ---------------------------------------------------------------------


      end subroutine calc_every_step
! *********************************************************************************


! *********************************************************************************
      subroutine no_radiation
! *********************************************************************************
 ! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('no_radiation')
      endif
! ---------------------------------------------------------------------

      rres(1:r%upper)=0.; rres(2) = 99. ! NO RADIATION

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('no_radiation')
      if (time_report) call write_clock_time ('no_radiation')
! ---------------------------------------------------------------------

      end subroutine no_radiation
! *********************************************************************************


! *********************************************************************************
      subroutine sloped_surfaces
! *********************************************************************************
 ! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('sloped_surfaces')
      endif
! ---------------------------------------------------------------------

! CALCULATE BEAM RADIATION, TOTAL RADIATION AND INCIDENCE ANGLE FOR EACH is-th PLANE
      axslp=r%slopes(ip)*rdconv
      axazm=r%axes(ip)*rdconv
      sinasl=sin(axslp)
      ip = ip+1
      cosasl=cos(axslp); cosasl=sign(amax1(abs(cosasl),1.e-06),cosasl)
      tanasl = sinasl/cosasl
      alf = sazm - axazm; alf = sign(amax1(abs(alf),1.e-06),alf)
! KEEP THE DIFFERENCE OF SOLAR AZIMUTH AND AXIS AZIMUTH BETWEEN 180 AND -180 DEGREES
      if(abs(alf)>pi) alf = alf - sign((2.*pi),alf)
      costtp = cosasl*coszen + sinasl*sinzen*cos(alf)
      costtp = sign(amax1(abs(costtp),1.e-06),costtp)

      tracking_mode: Select case (r%trackmode)
        case (1) ! FIXED SURFACE
          slope = axslp; azm = axazm
          
        case (2) ! SINGLE-AXIS TRACKER
          aux=abs(amod(axslp*dgconv,180.))
          if(aux>0.1.and.aux.le.179.9) then 
            slope = axslp; azm = sazm ! VERTICAL AXIS
          else
            azm = axazm + .5*sign(pi,alf)
            tanslp = tanzen*cos(sazm - azm); slope = atan(tanslp)
            if(slope<0.) slope = pi + slope
          endif 

        case (3) ! SLOPED AXIS

		!Horizontal axis		
		axazm_angle = r%axes(ip-1)	!Azimuth angle in degrees
		!East-west case
		if ((axslp == 0.d0) .and. ((axazm_angle  == 90.0) .or. (axazm_angle == - 90.0))) then
			slope = atan( tan(zenith)*abs( cos(sazm)) )
			if (abs(sazm) < 90.0) azm = 0.0
			if (abs(sazm) > 90.0) azm = pi
		!North-south case		
		elseif ((axslp == 0.d0) .and. ((axazm_angle  == 0.0) .or. (axazm_angle == pi))) then		
			if (sazm > 0.d0) azm = pi/2.0
			if (sazm < 0.d0) azm = -pi/2.0
			slope = atan( tan(zenith)*abs( cos(azm - sazm)) )	
		!Tilted axis			
		else
          tangam = sinzen*sin(alf)/sinasl/costtp
          gam = atan(tangam); 
		  if(gam<.0.and.alf>0.) gam = pi+gam
          if(gam>.0.and.alf<0.) gam = gam-pi
          tanslp = tanasl/cos(gam); slope = atan(tanslp)
          if(slope<0.) slope = pi+slope; azm = gam + axazm
          if(abs(azm)>pi) azm = azm - sign((2.*pi),azm)
          endif 

        case (4) ! TWO-AXIS TRACKER
          azm = sazm; slope = zenith
      end select tracking_mode    

      cosslp = cos(slope); sinslp = sin(slope)
      costt = min(1.,cosslp*coszen+sinslp*sinzen*cos(sazm-azm))
      theta = acos(costt)*dgconv; slope = slope*dgconv

! BEAM AND GROUND REFLECTED RADIATION INDEPENDENT OF TILTED SURFACE MODEL
      rb = amax1(costt,0.)/coszen; hbeam = hb*rb
      hgrf  = hhor*r%rho*.5*(1.-cosslp)
     
      tilted_surface_skymode: select case (r%skymode)
        case (1)               !ISOTROPIC SKY MODEL FOR TILTED SURFACE DIFFUSE 
          hdiff = hd*.5*(1.+cosslp)

        case (2)                         !HAY MODEL FOR TILTED SURFACE DIFFUSE
          ai = hb/hextra; hdiff = hd*(.5*(1.-ai)*(1.+cosslp)+ai*rb)

        case (3)                                  !REINDL TILTED SURFACE MODEL
          ai = hb/hextra
          f  = sqrt(hb/max(hhor,1.0)); scube = (sin(slope*0.5*rdconv))**3
          hdiff = hd*(0.5*(1.-ai)*(1.+cosslp)*(1.+f*scube)+ai*rb)

        case (4)           !PEREZ POINT SOURCE MODEL (SANDIA REPORT OCT, 1988)
          hdn = hb/coszen
          if ( hd < 0.00001 ) then
            eps = 99999.
          else
            epsiln = ( hd + hdn ) / hd
            eps = ( epsiln+1.041*zenith**3 )/( 1.+1.041*zenith**3 )
          endif
          skyb = hd/hextra
          if (eps>0.0 ) nbin = 1; if (eps>1.065) nbin = 2
          if (eps>1.23) nbin = 3; if (eps>1.5  ) nbin = 4
          if (eps>1.95) nbin = 5; if (eps>2.8  ) nbin = 6
          if (eps>4.5 ) nbin = 7; if (eps>6.2  ) nbin = 8

          p1=p11(nbin)+p12(nbin)*skyb+p13(nbin)*zenith;p1=max(p1,.0)
          p2=p21(nbin)+p22(nbin)*skyb+p23(nbin)*zenith
          a1=max(costt,0.); b1=amax1(cos(85.0*rdconv),coszen)

          if (hb > 0.) then
            hdiff = hd*(0.5*(1.-p1)*(1.+cosslp)+p1*a1/b1+p2*sinslp)
          else        !DIFFUSE ASSUMED ISOTROPIC IF BEAM NOT POSITIVE
            hdiff = hd*0.5*(1.+cosslp)
          endif 
          hdiff = max(hdiff,0.)

      end select tilted_surface_skymode

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('sloped_surfaces')
      if (time_report) call write_clock_time ('sloped_surfaces')
! ---------------------------------------------------------------------

      end subroutine sloped_surfaces
! *********************************************************************************

! ***************************************************************************
      subroutine suncalc
! ***************************************************************************
! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) then
        call get_clock_on_entry('suncalc')
      endif
! ---------------------------------------------------------------------

!CALCULATE SUN POSITION      
      sazm = 0. 
      coszen=cc*coshr+ss; coszen=sign(amax1(abs(coszen),1.e-06),coszen)
      zenith=acos(coszen); sinzen=sin(zenith); tanzen=sinzen/coszen

      if (abs(sinzen) .ge. 1e-06) then
        sinazm=cosdec*sinhr/sinzen
        sinazm=sign(amin1(abs(sinazm),1.),sinazm); sazm = asin(sinazm)
       !DETERMINE WETHER THE SOLAR AZIMUTH IS GREATER THAN 90 DEGREES BY 
       !COMPARING THE HOUR ANGLE WITH THE HOUR ANGLE AT WHICH THE SOLAR 
       !AZIMUTH IS +/- 90 DEGREES
        cwew = tandec/r%tanlat; cwew=sign(amin1(abs(cwew),1.),cwew)
        wew = pi; if(r%latit*(decl-r%latit).le.0.) wew = acos(cwew)
        if((abs(w)-abs(wew))*r%latit*(decl-r%latit).le.0.) sazm = sign(pi,sazm) - sazm
       !DON'T ALLOW THE SOLAR AZIMUTH TO BE GREATER THAN 180 DEGREES.
        if(abs(sazm) > pi) sazm = sazm - sign((2.*pi),sazm)
      endif

! type109_time_report: 24.02.2005 TW ------------------------------
      if (time_report) call get_clock_on_exit('suncalc')
      if (time_report) call write_clock_time ('suncalc')
! ---------------------------------------------------------------------

      end subroutine suncalc
! ***************************************************************************


      end subroutine rad_proc
! ***************************************************************************