ALGORITHMIC_MODEL Schedule ABSTRACT "Fetch data from SCHEDULE. Handle leap yrs & DST. Generate events. This is the algorithmic schedule object for IDA CE. Leap years and DST are handled correctly. Discontinuities are handled, both at midnight, and when jumps occur with double time coordinates. Smoothing of each 24hr profile can be done in the parameter processing. It is controlled by two parameters sample_width Width of bins that are redistributed by a symmetric staircase function to adjacent bins. smoothing_fanout Number of bins on each side that are affected. Example: When fanout is 2, a bin with function value 1 is spread as | | 3/9 --- Sum under stair is | | | | 2/9 --- --- 9 * 1/9 = 1 | | | | 1/9 --- --- | | | ----------------------------------------------- | | original bin location Smoothing_fanout = 0 suppresses smoothing. Daylight saving During a yearly simulation, integration time will follow standard clock time. In the daylight saving period, this time is increased by one hour (in UT2Calendar), when schedules are accessed. If the integration starts within the dayligth saving period (and extends over one or more full days), the integration time should start at 23 hrs, standard clock time. ! Planned modification: Specify DST start and stop days in one of three ways: > 0 Use days as specified, for all years 0 Calculate according to EU standard (last sundays in Mar-Oct), depending on each current year -1 No daylight saving _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ Revisions: LE 020305 Removed bug in case where time step reaches over PG midnight. Event before midnight was not detected. PG 020301 Function EventTime now requires true time from the beginning of the year LE 020220 Function StepF moved to ut2calendar.nmf LE 011218 Bug in Table_Lookup PG 011026 OutSignal renamed to OutSignalLink PG 010330 <29> Day number for rules: cal_day = mmdd MV 010322 New link type is added LE 990224 Bug in Table_Lookup fixed AB 981207 Time event bug fixed (short pulses were missed) AB 980514 DST on & off fixed AB 980510 Midnight event rectified AB 980121 DST_hours corrected for dayno specified as 1.. AB 980114 Error if not calendar time AB 971203 Documentation updated MV 971125 Smoothing_Fanout default is changed to zero AB 971020 Revised ut2calendar; removed unused parameters. LE 971016 Bug fixes in ut2calendar, added 'days' = day no, from 1900-01-01 = 1 AB 970928 'Time' in 'universal time', i.e. seconds from 1900-01-01. AB 970927 Integration time is in hours from beginning of 1st year, but 'Time' is always delivered in seconds. AB 970831 Added smoothing AB 970821 Int replaces Nint in days:= day_in_year 1-365.5 PG 970624, 970820 HH 970506 " EQUATIONS /* Time is the built in time variable, always in seconds; 0 at beginning of year 1900. */ /* Seconds from the beginning of year 1900 */ seconds := Time; /* Enter universal time seconds [s], get calendar data back */ CALL UT2Calendar(seconds, DST_hour_1, DST_hour_2, year_no, day_in_year, hour1_in_year, hour1, month, day, weekday); If year_no <= 0 THEN CALL nmf_error ("Integration time must use calendar format"); END_IF; /* The valid rule */ /* <29> = rule_no := Hit_Rule(n, ndpw, year_no, day_in_year, weekday, Case_Date, */ rule_no := Hit_Rule(n, ndpw, year_no, 100*month+day, weekday, Case_Date, Start_Date, End_Date, Wd_Code); /* If new rule, restart table-lookup */ IF rule_no>0 THEN tab_ptr := IF rule_no != prev_rule_no THEN tab_beg[rule_no] ELSE prev_tab_ptr END_IF; END_IF; /* The output variable */ IF hour1 < prev_hour THEN /* Check if midnight passed */ /* lookup should be done, also in end of previous day to check for time events */ CALL Table_Lookup(nTab,tab_time,tab_value,prev_tab_ptr,24.0, val_24,disc_time); IF disc_time >= 0.0 THEN CALL EventTime (1.0, DST_Hour((day_in_year-1)*24+disc_time, DST_hour_1, DST_hour_2)*3600) ELSE val_0 := IF rule_no>0 THEN value_col[1, rule_no] ELSE 0.0 END_IF; IF abs (val_0 - val_24) > 1e-6 THEN CALL EventTime (1.0, DST_Hour(day_in_year*24, DST_hour_1, DST_hour_2)*3600) END_IF; END_IF; END_IF; CALL Table_Lookup(nTab,tab_time,tab_value,tab_ptr,hour1, curr_val,disc_time); Out_Var := IF rule_no>0 THEN In_Par[rule_no] * curr_val ELSE 0.0 END_IF; /* Suppress scheduled events in hour skipped when DST goes on */ IF disc_time > 0 AND (prev_hr1_in_yr > DST_hour_1 OR hour1_in_year < DST_hour_1) THEN /* = CALL EventTime (1.0, (day_in_year*24+disc_time)*3600) */ CALL EventTime (1.0, DST_Hour(day_in_year*24+disc_time, DST_hour_1, DST_hour_2)*3600) END_IF; /* Generate events when DST goes on and off */ /* Jump event at DST on */ IF prev_hr1_in_yr < DST_hour_1 AND hour1_in_year > DST_hour_1+1 THEN /* = CALL EventTime (1.0, (DST_hour_1+1)*3600) */ CALL EventTime (1.0, DST_hour_1*3600) END_IF; /* Knee events around two hours with half pace, prior to DST off */ IF prev_hr1_in_yr < DST_hour_2-1 AND hour1_in_year > DST_hour_2-1 THEN /* = CALL EventTime (2.0, (DST_hour_2-1)*3600) */ CALL EventTime (2.0, (DST_hour_2-2)*3600) END_IF; IF prev_hr1_in_yr < DST_hour_2 AND hour1_in_year > DST_hour_2 THEN CALL EventTime (2.0, DST_hour_2*3600) END_IF; /* Update state memory */ prev_rule_no := rule_no; prev_hr1_in_yr := hour1_in_year; prev_hour := hour1; prev_tab_ptr := tab_ptr; LINKS /* type name variables... */ Generic Out_Link Out_Var; Generic OutSignalLink Out_Var; VARIABLES /* type name role [def [min max]] description */ Generic Out_Var OUT "output variable" Generic seconds LOC "second counter" Generic disc_time LOC "time of possible event [hr]" Generic weekday LOC "day in week counter" Generic year_no LOC "years from start year" Generic month LOC "month in current year (1-12)" Generic day LOC "day in current month (1-31)" Generic hour1_in_year LOC "hour in current year" Generic prev_hr1_in_yr A_S 9999 0 BIG "previous ditto" Generic hour1 LOC "hour in current day (0-24)" Generic day_in_year LOC "day in current year" Generic rule_no LOC "the number of the valid rule" Generic prev_Rule_no A_S 0 "preious rule number" Generic prev_hour A_S "previous time coord" factor curr_val LOC "value from table lookup" factor val_0 LOC "1st value in current rule" factor val_24 LOC "last value in previous rule" factor tab_ptr LOC "pointer into table" factor prev_tab_ptr A_S 1 1 BIG "previous table pointer" MODEL_PARAMETERS /* type name role [def min max] description */ INT ndpw SMP 7 7 7 "number of days per week" INT n SMP 1 1 BIGINT "number of rules" INT len SMP 1 1 BIGINT "length of a profile" INT nTab CMP 2 2 BIGINT "no of table cells" INT nS CMP 2 2 BIGINT "no of sample bin" PARAMETERS /* type name role [def min max] description */ /* Generic Second_1 S_P 0 0 31536000 "start second for the simulation in the year" Generic Year_1 S_P 1997 0 BIG "start year for the simulation" */ Generic DST_start S_P 88 0 366 "start day for daylight savings time" Generic DST_end S_P 298 0 366 "end day for daylight savins time" Generic DST_hour_1 C_P "start hour for DST" Generic DST_hour_2 C_P "end hour for DST" Factor Sample_Width S_P 0.2 SMALL BIG "smoothing resolution, width of sample bin" Factor Smoothing_Fanout S_P 0 0 BIG "no of samples on each side, affected by sample" Generic In_Par[n] S_P 1.0 "values for rules" Generic Wd_Code[ndpw, n] S_P "weekday code for rule i" Generic Start_Year[n] S_P -1 "start year for rule i" Generic End_Year[n] S_P -1 "end year for rule i" Generic Start_Month[n] S_P -1 -1 12 "start month for rule i" Generic End_Month[n] S_P -1 -1 12 "end month for rule i" Generic Start_Day[n] S_P -1 -1 31 "start day for rule i" Generic End_Day[n] S_P -1 -1 31 "end day for rule i" Factor Tab_Beg[n] C_P 1 1 BIG "beg index in table" Generic time_col[len, n] S_P "time cols for profiles" Generic value_col[len, n] S_P "value cols for profiles" Generic Case_Date[n] C_P "code of restrictions: 1 - no restrictions; 2 - from ... to every year 3 - .. to and from .. every year 4 - from ... to with fixed years" Generic Start_Date[n] C_P "start day in year for rule i" Generic End_Date[n] C_P "end day in year for rule i" /* smoothed profiles, stored after each other in combined vectors */ Factor Tab_Time[nTab] C_P "time vector" Factor Tab_Value[nTab] C_P "value vector" Factor bin_area[nS] C_P "areas of sample bins" PARAMETER_PROCESSING /* size of bin vector */ nS := Nint (24.0/sample_width + 1); /* extend max profile to infinities */ nTab := n * (nS + 2); DST_hour_1 := 24.0 * (DST_start - 1) + 2.0; DST_hour_2 := 24.0 * (DST_end - 1) + 3.0; CALL Smooth_Rules(len, n, time_col,value_col, smoothing_fanout, nTab, tab_time,tab_value, tab_beg, nS, bin_area); CALL Prepare_Rules(n, Start_Year, End_Year, Start_Month, End_Month, Start_Day, End_Day, Case_Date, Start_Date, End_Date) /* ======== The internal functions =========== */ /* ********************************************************************* CALL Smooth_Rules(len, n, time_col,value_col, smoothing_fanout, nTab, tab_time,tab_value, tab_beg, nS, bin_area); Smooth_Rules Input len,n : size of supplied rule matrices. time_col,value_col : supplied rule matrices, one column for each rule. smoothing_fanout : no of bins on each side affected by sample, each sample staple is replaced by a staircase with width = (2*smoothing_fanout+1)*sample_width. nTab : length of smoothed table vector. Output tab_time,tab_value : smoothed rules concatenated. tab_beg : pointers to rules within table. */ FUNCTION VOID Smooth_Rules(len, n, time_col,value_col, smoothing_fanout, nTab, tab_time,tab_value, tab_beg, nS, area) LANGUAGE F77 INPUT INT len, n, nTab, nS; FLOAT time_col[len, n], value_col[len, n], smoothing_fanout; OUTPUT FLOAT tab_time[nTab], tab_value[nTab], tab_beg[n], area[nS]; CODE SUBROUTINE Smooth_Rules(len, n, time_col,value_col, & smoothing_fanout, & nTab, tab_time,tab_value, tab_beg, & nS, area) INTEGER len, n, nS, nTab DOUBLE PRECISION time_col(1:len, 1:n),value_col(1:len, 1:n), & smoothing_fanout, & tab_time(nTab), tab_value(nTab), tab_beg(n), & area(0:nS-1) INTEGER fanout, k, row, rule, s, tab_i REAL*8 inf, & sample_area, tot_area, & neighbor_lot, unit_lot, & slope, & t1_a,t2_a, t1_p,t2_p, t1_s,t2_s, ! a:subarea p:profile s:sample & v1_a,v2_a, v1_p,v2_p, v1_s,v2_s ctest DATA inf /1E30/ DATA inf /25/ fanout = Nint(smoothing_fanout) tab_i = 1 * Suppress smoothing if smoothing_fanout is zero IF (fanout .EQ. 0) GOTO 2000 * Scan rules in 'time_col,value_col', for each rule, * add smoothed version to end of 'tab_time,tab_value' DO 1000 rule = 1, n * Clear sample bins DO 100 s = 0, nS-1 area(s) = 0 100 CONTINUE * Scan samples, calculate profile area for each sample bin t1_p = -inf v2_p = value_col(1, rule) t2_p = time_col(1, rule) row = 1 tot_area = 0 DO 200 s = 0, nS-1 * set limits of current bin t1_s = 24.0 * (s - 0.5)/(nS-1) t2_s = 24.0 * (s + 0.5)/(nS-1) IF (s .EQ. 0) THEN t1_s = 0 ELSE IF (s .EQ. nS-1) THEN t2_s = 24 ENDIF sample_area = 0 t1_a = t1_s * repeat from here when profile segment ends inside sample bin * scan profile to locate current subarea 120 CONTINUE IF (t1_a .GE. t2_p) THEN t1_p = t2_p v1_p = v2_p row = row + 1 t2_p= time_col(row, rule) v2_p = value_col(row, rule) GOTO 120 ENDIF * now, t1_a is inside profile segment * t1_p <= t1_a < t2_p slope = (v2_p - v1_p) / (t2_p - t1_p) v1_a = v1_p + (t1_a - t1_p) * slope IF (t2_p .LT. t2_s) THEN t2_a = t2_p v2_a = v2_p ELSE t2_a = t2_s v2_a = v1_p + (t2_a - t1_p) * slope ENDIF sample_area = sample_area + & (v1_a + v2_a)/2 * (t2_a - t1_a) * check if right end of sample bin has been reached t1_a = t2_a IF (t2_s .GT. t2_a) GOTO 120 tot_area = tot_area + abs (sample_area) * spread area onto adjacent bins unit_lot = sample_area / (24.0/(nS-1)) / (fanout+1)**2 DO 160 k = -fanout, fanout neighbor_lot = unit_lot * (fanout+1-iabs(k)) IF (s+k .LT. 0 .OR. s+k .GT. nS-1) THEN area(s) = area(s) + neighbor_lot ELSE area(s+k) = area(s+k) + neighbor_lot ENDIF 160 CONTINUE 200 CONTINUE *. end samples * Compress profile, telescoping linear stretches. * t1_a,v1_a store beginning of linear stretch. * Slope is valid for current linear strech, and is updated for each * point that extends the stretch tab_beg(rule) = tab_i t1_a = -inf v1_a = area(0) t2_s = 0 v2_s = v1_a slope = 0 DO 900 s = 0, nS-1 * fetch profile value (v2_s) for new sample (at t2_s), * calculate corresponding value (v2_a) * as extrapolated along current slope. t1_s = t2_s v1_s = v2_s t2_s = 24.0 * s / (nS - 1) v2_s = area(s) v2_a = v1_a + slope * (t2_s - t1_a) IF (abs (v2_a - v2_s) .GT. 1E-4 * tot_area & .OR. s .EQ. nS-1) THEN * store previous stretch, prepare for new tab_time(tab_i) = t1_a tab_value(tab_i) = v1_a tab_i = tab_i + 1 t1_a = t1_s v1_a = v1_s ENDIF * adjust slope slope = (v2_s - v1_a) / (t2_s - t1_a) 900 CONTINUE tab_time(tab_i) = inf tab_value(tab_i) = v1_a tab_i = tab_i + 1 1000 CONTINUE RETURN * Move profiles to common vector without smoothing, * when smoothing_fanout = 0 2000 CONTINUE DO 2200 rule = 1, n tab_beg(rule) = tab_i t2_p = -inf tab_time(tab_i) = t2_p v2_p = value_col(1, rule) tab_value(tab_i) = v2_p tab_i = tab_i + 1 DO 2100 row = 1, len t1_p = t2_p v1_p = v2_p t2_p = time_col(row, rule) v2_p = value_col(row, rule) IF (t2_p .LE. t1_p & .AND. t2_p .GE. 24) GOTO 2100 tab_time(tab_i) = t2_p tab_value(tab_i) = v2_p tab_i = tab_i + 1 2100 CONTINUE tab_time(tab_i) = inf tab_value(tab_i) = v1_p tab_i = tab_i + 1 2200 CONTINUE RETURN END END_CODE END_FUNCTION /* ********************************************************************* Table lookup Input tab_time : Time vector; monotonously growing; 1st value less than any expected input, last value greater than any expected input. tab_value : Value vector; parallel with tab_time. tab_ptr : start index in vectors; could point after current time. arg : time of day [hr] Output curr_val : value at arg tab_ptr : updated to point at either end of relevant interval. disc_time : time of detected time event [hr]; negative if no event detected. Example tab_time -inf 0 8 9 11 13 15 17 18 24 inf tab_value 0 0 0 20 22 21 14 11 0 0 0 Active table may be part of longer vector, with several tables concatenated. */ FUNCTION VOID Table_Lookup(nTab,tab_time,tab_value,tab_ptr, arg, curr_val,disc_time) LANGUAGE F77 INPUT INT nTab; FLOAT tab_time[nTab], tab_value[nTab], tab_ptr, arg; OUTPUT FLOAT curr_val, tab_ptr, disc_time; CODE SUBROUTINE Table_Lookup ( & nTab,tab_time,tab_value,tab_ptr, arg, & curr_val, disc_time) INTEGER nTab DOUBLE PRECISION tab_time(nTab), tab_value(nTab), tab_ptr, arg, & curr_val, disc_time INTEGER iTab_ptr DOUBLE PRECISION val, slope * INTERPOLATE LINEARLY IN TABLE 'tab_Time,tab_Value', * find values 'val' & 'slope' corresponding to argument 'arg'. * 'tab_ptr' points into table and is adjusted, * up or down, until a suitable interval is found. * * 'disc_time' > 0 points to earliest event found while adjusting * 'tab_ptr', < 0 if no event. * * It is assumed that the outmost segments stretch to +- 'infinity'. * 'arg' is asserted to stay inside this interval. * REAL*8 inf, t1,t2, v1,v2 ctest DATA inf /1E30/ DATA inf /25/ * get integer pointer iTab_ptr = nint (tab_ptr) IF (arg .LE. -inf .OR. arg .GE. inf) THEN CALL nmf_error ('Invalid time argument in Table_Lookup') ENDIF t1 = tab_time(iTab_ptr) C9811 disc_time = -1.0 disc_time = inf IF (arg .LT. t1) THEN * search backwards; t2 < t1 20 CONTINUE IF (iTab_ptr .GT. 1) THEN iTab_ptr = iTab_ptr - 1 t2 = tab_time(iTab_ptr) IF (t2 .GT. t1) THEN CALL nmf_error ('Poorly organized table in Table_Lookup') C9811 ELSE IF (t2 .GE. t1 .AND. t2 .LT. disc_time) THEN disc_time = t2 ENDIF IF (arg .LT. t2) THEN t1 = t2 GOTO 20 ENDIF C9902 LE- ELSE t2 = t1 t1 = tab_time(2) IF (t2 .GT. t1) THEN CALL nmf_error ('Poorly organized table in Table_Lookup') ELSE IF (t2 .GE. t1 .AND. t2 .LT. disc_time) THEN disc_time = t2 ENDIF C9902 LE. ENDIF v1 = tab_value(iTab_ptr+1) ELSE * search forwards; t2 > t1 40 IF (iTab_ptr .LT. nTab) THEN iTab_ptr = iTab_ptr + 1 t2 = tab_time(iTab_ptr) IF (t2 .LT. t1) THEN CALL nmf_error ('Poorly organized table in Table_Lookup') C9811 ELSE IF (t2 .LE. t1) THEN ELSE IF (t2 .LE. t1 .AND. t2 .LT. disc_time) THEN disc_time = t2 ENDIF IF (arg .GT. t2) THEN t1 = t2 GOTO 40 ENDIF C9902 LE- ELSE t2 = t1 t1 = tab_time(nTab-1) IF (t2 .LT. t1) THEN CALL nmf_error ('Poorly organized table in Table_Lookup') ELSE IF (t2 .LE. t1 .AND. t2 .LT. disc_time) THEN disc_time = t2 ENDIF C9902 LE. ENDIF v1 = tab_value(iTab_ptr-1) ENDIF v2 = tab_value(iTab_ptr) C0112 LE+ if (t1 .eq. t2) then val = v2 c write(*,*) 'Table_lookup: arg ', arg, disc_time c write(*,*) 'nTab, iTab_ptr ',nTab, iTab_ptr c write(*,*) 't1, t2 ', t1, t2 c write(*,*) 'v1, v2 ', v1, v2 else C0112 LE. slope = (v2 - v1) / (t2 - t1) val = (arg - t1) * slope + v1 endif curr_val = val Tab_ptr = iTab_ptr C9811 C9902 IF (disc_time .EQ. inf) disc_time = -1 IF (disc_time .EQ. inf .OR. disc_time .GT. arg) disc_time = -1 RETURN END END_CODE END_FUNCTION /* ********************************************************************* Convert the date restrictions to easy-to-use format. START_DATE & END_DATE are day numbers in the year (1-365.5). Conversion is done without regard to leap years. */ FUNCTION VOID Prepare_Rules(n, Start_Year, End_Year, Start_Month, End_Month, Start_Day, End_Day, Case_Date, Start_Date, End_Date) LANGUAGE F77 INPUT INT n; FLOAT Start_Year[n], End_Year[n], Start_Month[n], End_Month[n], Start_Day[n], End_Day[n]; OUTPUT FLOAT Case_Date[n], Start_Date[n], End_Date[n]; CODE CAB SUBROUTINE PREPARE_RULES( N, YEAR, START_YEAR, END_YEAR, SUBROUTINE PREPARE_RULES( N, START_YEAR, END_YEAR, & START_MONTH, END_MONTH, START_DAY, END_DAY, & CASE_DATE, START_DATE, END_DATE) INTEGER N,J CAB DOUBLE PRECISION YEAR, TMP_YEAR DOUBLE PRECISION TMP_YEAR DOUBLE PRECISION START_YEAR(N), END_YEAR(N), START_MONTH(N) DOUBLE PRECISION END_MONTH(N), START_DAY(N), END_DAY(N) DOUBLE PRECISION CASE_DATE(N), START_DATE(N), END_DATE(N) DOUBLE PRECISION SHIFT_DAYS(12), LAST_DAYS(12) C<29> - C DATA SHIFT_DAYS /0, 31, 59.5, 90.5, 120.5, 151.5, 181.5, C & 212.5, 243.5, 273.5, 304.5, 334.5 / C DATA LAST_DAYS /31, 28.5, 31, 30, 31, 30, 31, C & 31, 30, 31, 30, 31 / C<29> . C SAVE SHIFT_DAYS DO 30 J=1,N IF (START_MONTH(J) .LE. 0.0) START_MONTH(J) = 1 IF (START_DAY(J) .LE. 0.0) START_DAY(J) = 1 IF (END_MONTH(J) .LE. 0.0) END_MONTH(J) = 12 C<29> - C IF (END_DAY(J) .LE. 0.0) C & END_DAY(J) = LAST_DAYS(NINT(END_MONTH(J))) C START_DATE(J) = START_DAY(J) + SHIFT_DAYS(NINT(START_MONTH(J))) C END_DATE(J) = END_DAY(J) + SHIFT_DAYS(NINT(END_MONTH(J))) C<29> + IF (END_DAY(J) .LE. 0.0) END_DAY(J) = 50 START_DATE(J) = START_DAY(J) + 100*START_MONTH(J) - 0.1 END_DATE(J) = END_DAY(J) + 100*END_MONTH(J) + 0.1 C<29> . IF (START_YEAR(J) .LE. 0.0) THEN IF (END_YEAR(J) .LE. 0.0) THEN * No year specified C<29> - C IF (START_DATE(J) .LT. 1.5 .AND. C & END_DATE(J) .GE. 365.0) THEN C<29> + IF (START_DATE(J) .LE. 101 .AND. & END_DATE(J) .GE. 1231) THEN C<29> + CASE_DATE(J) = 1.0 ELSEIF (START_DATE(J) .LE. END_DATE(J)) THEN CASE_DATE(J) = 2.0 ELSE CASE_DATE(J) = 3.0 ENDIF ELSE * Only end year specified, set start year reative to end year IF (START_DATE(J) .LT. END_DATE(J)) THEN TMP_YEAR = END_YEAR(J) ELSE TMP_YEAR = END_YEAR(J) - 1.0 ENDIF CAB START_DATE(J) = 1000.0 * (TMP_YEAR-YEAR) + START_DATE(J) CAB END_DATE(J) = 1000.0 * (END_YEAR(J)-YEAR) + END_DATE(J) C<29>: 1000 -> 10000 START_DATE(J) = 10000.0 * (TMP_YEAR-1900) + START_DATE(J) END_DATE(J) = 10000.0 * (END_YEAR(J)-1900) + END_DATE(J) CASE_DATE(J) = 4.0 ENDIF ELSE IF (END_YEAR(J) .LE. 0.0) THEN * Only start year specified, set end year reative to start year IF (START_DATE(J) .LT. END_DATE(J)) THEN TMP_YEAR = START_YEAR(J) ELSE TMP_YEAR = START_YEAR(J) + 1.0 ENDIF START_DATE(J) = CAB & 1000.0 * (START_YEAR(J)-YEAR) + START_DATE(J) CAB END_DATE(J) = 1000.0 * (TMP_YEAR-YEAR) + END_DATE(J) & 10000.0 * (START_YEAR(J)-1900) + START_DATE(J) END_DATE(J) = 10000.0 * (TMP_YEAR-1900) + END_DATE(J) ELSE * Both start and end year specified START_DATE(J) = CAB & 1000.0 * (START_YEAR(J)-YEAR) + START_DATE(J) CAB END_DATE(J) = 1000.0 * (END_YEAR(J)-YEAR) + END_DATE(J) & 10000.0 * (START_YEAR(J)-1900) + START_DATE(J) END_DATE(J) = 10000.0 * (END_YEAR(J)-1900) + END_DATE(J) ENDIF CASE_DATE(J) = 4.0 ENDIF 30 CONTINUE RETURN END END_CODE END_FUNCTION /* ********************************************************************* A function for finding a valid rule number <29> = date = day in year (1-366) date = day in year = mmdd */ FUNCTION FLOAT Hit_Rule(n, nwd, rel_year, date, weekday, Case_Date, Start_Date, End_Date, Wd_Code) LANGUAGE F77 INPUT INT n, nwd; FLOAT Case_Date[n], Start_Date[n], End_Date[n], Wd_Code[nwd,n], rel_year, date, weekday; CODE FUNCTION HIT_RULE (N, NWD, REL_YEAR, DATE, WEEKDAY, & CASE_DATE, START_DATE, END_DATE, WD_CODE) INTEGER N,NWD,I DOUBLE PRECISION HIT_RULE, REL_YEAR, DATE, WEEKDAY DOUBLE PRECISION TDATE, WD_CODE(NWD, N) DOUBLE PRECISION CASE_DATE(N), START_DATE(N), END_DATE(N) HIT_RULE = 0.0 DO 100 I=1,N IF ( WD_CODE( NINT(WEEKDAY), I) .EQ. 0.0) GOTO 100 GOTO (10,20,30,40) NINT( CASE_DATE(I)) 10 HIT_RULE = I GOTO 200 20 IF (DATE .GE. START_DATE(I) .AND. & DATE .LE. END_DATE(I)) GOTO 10 GOTO 100 30 IF (DATE .LE. START_DATE(I) .OR. & DATE .GE. END_DATE(I)) GOTO 10 GOTO 100 40 TDATE = 10000.0 * REL_YEAR + DATE IF (TDATE .GE. START_DATE(I) .AND. & TDATE .LE. END_DATE(I)) GOTO 10 100 CONTINUE 200 CONTINUE RETURN END END_CODE END_FUNCTION /* ********************************************************************* A function for evaluating the profile using linear interpolation. The data columns are assumed to be sorted in rising order. Discontinuities are allowed but not handled. */ FUNCTION FLOAT Eval_Prof(len, n, j, col_1, col_2, time_24, values) LANGUAGE F77 INPUT INT len, n, j; FLOAT col_1[ len, n], col_2[ len, n], time_24, values[n]; CODE FUNCTION EVAL_PROF( LEN, N, J, COL_1, COL_2, TIME_24, VALUES) INTEGER LEN, N, I, J DOUBLE PRECISION EVAL_PROF, COL_1, COL_2, TIME_24, VALUES DOUBLE PRECISION DELTA_T, DELTA_V, SLOPE, RESULT DIMENSION COL_1(1:LEN, 1:N), COL_2(1:LEN, 1:N), VALUES(1:N) RESULT = 0.0 IF (TIME_24 .LE. COL_1( 1, J)) THEN RESULT = COL_2( 1, J) GOTO 20 ENDIF DO 10 I=2,LEN IF (TIME_24 .LT. COL_1( I, J)) THEN DELTA_T = COL_1( I, J) - COL_1( I-1, J) DELTA_V = COL_2( I, J) - COL_2( I-1, J) IF ( ABS( DELTA_T) .LT. 0.0000001) GOTO 10 SLOPE = DELTA_V / DELTA_T RESULT = COL_2( I-1, J) + (TIME_24 - COL_1( I-1, J)) * SLOPE GOTO 20 ENDIF 10 CONTINUE RESULT = COL_2( LEN, J) 20 CONTINUE EVAL_PROF = RESULT * VALUES(J) RETURN END END_CODE END_FUNCTION /* ********************************************************************* A function for checking leap years. */ FUNCTION FLOAT Leap_year( year) LANGUAGE F77 INPUT FLOAT year; CODE FUNCTION LEAP_YEAR( YEAR) DOUBLE PRECISION LEAP_YEAR, YEAR, RESULT RESULT = 0.0 IF (MOD( NINT( YEAR), 4) .EQ. 0) THEN RESULT = 1.0 ENDIF IF (MOD( NINT( YEAR), 100) .EQ. 0) THEN RESULT = 0.0 ENDIF IF (MOD( NINT( YEAR), 400) .EQ. 0) THEN RESULT = 1.0 ENDIF LEAP_YEAR = RESULT RETURN END END_CODE END_FUNCTION /* ********************************************************************* A system function to signal time dependent events. */ FUNCTION VOID EventTime(iventType, iventTime) LANGUAGE F77 INPUT FLOAT iventType, iventTime; EXTERNAL EventTime (iventType, iventTime) END_FUNCTION END_MODEL