[TRNSYS-users] Type109 high latitude bug
Janne Paavilainen
jip at du.se
Mon Jun 25 08:47:03 PDT 2007
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 at 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 at byg.dtu.dk <mailto:jif at byg.dtu.dk>
______________________________________
-------------- next part --------------
! 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
! ***************************************************************************
More information about the TRNSYS-users
mailing list