! 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: 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)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.'') then call MESSAGES(217,'','fatal',INFO(1),INFO(2)) return endif do while (part /= '') !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 select case (part) !read definitions case ('') 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 ('') 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 ('') 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 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)0 .and. ifile(3) 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)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. xtime0) 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)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)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 (nlastdayu(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 (suns0) 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=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 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 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 (uk3) 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.w1pi) 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 ! ***************************************************************************