Определить количество дней между двумя заданными датами. Фортран 77

C          Calculate conforming to complex calendarical contortions.
C   Astronomer Simon Newcomb determined that the tropical year of 1900
c contained 31556925.9747 seconds, or 365.24219879 days.
c Subsequent definitions involve "no measurable differences",
c whereas in 45BC (when the Julian calendar was adopted), the year was
c 365.24232 days long, going by modern calculations.
c   The "tropical" year is the time between the same equinoxes, and thus
c contains the effect of the precession of the Earth's axis, which would
c otherwise cause the seasons to likewise precess around the "fixed star"
c year, that being the time for midnight to point in the same direction
c amongst the "fixed" stars. Specifically, the vernal equinox
c (the northern hemisphere's spring equinox: cultural colonialism)
c is meant to hover around the 21'st of March, although it may fall
c within the 19'th or 20'th, which last has been the most popular in
c the 20'th century, until the leap year of 2000 resynchronised
c the civil calendar.
C   By contrast, the computations of astrologers are still based on the
c constellations as oriented in Babylonian times...
 
c   So, to add .24219879 to 365 days...
c   Adjustment per year   Nett     Discrepancy remaining.
c   +1/4     +.25        +.25      -.00780121  5h 48m 45.98s
c   -1/100   -.01         .24      +.00219879    -11m 14.02s
c   +1/400   +.0025       .2425    -.00030121        -26.02s
c   -1/4000  -.00025      .24225   -.00005121         -4.42s
c
c   The remnant of -.00005121 (meaning that the calendar year is too long)
c amounts to needing to drop one day on 19,527 years, and while this could be
c accommodated nicely enough by -1/20000 to give a calendar year of
c 365.24420 days with a remaining discrepancy of -.00000121 or -.01sec/year,
c there is a problem. The Earth's spin is slowing by a similar amount.
 
       TYPE DateBag     !Pack three parts into one.
        INTEGER YEAR,MONTH,DAY  !The usual suspects.
       END TYPE DateBag     !Simple enough.
 
       CHARACTER*9 MONTHNAME(12),DAYNAME(0:6)   !Re-interpretations.
       DATA MONTHNAME/"January","February","March","April","May","June",
     1  "July","August","September","October","November","December"/
       DATA DAYNAME /"Sunday","Monday","Tuesday","Wednesday","Thursday",
     1  "Friday","Saturday"/    !Index this array with DayNum mod 7.
 
       INTEGER*4 JDAYSHIFT  !INTEGER*2 just isn't enough.
       PARAMETER (JDAYSHIFT = 2415020)  !Thus shall 31/12/1899 give 0, a Sunday.
       DOUBLE PRECISION DAYSINYEAR  !A real ache.
       PARAMETER (DAYSINYEAR =  365.24219879D0) !The "D0" demands DOUBLE PRECISION precision.
       INTEGER*4 SECONDSINDAY   !This has its uses.
       PARAMETER (SECONDSINDAY = 24*60*60)  !86400. Disregarding "leap" seconds.
 
       PARAMETER NZBASEPLACE = "Mt. Cook Trig, Wellington." !Name the place.
       DOUBLE PRECISION NZBASELAT,NZBASELONG    !Alas, DATA statements do not allow arithmetic expressions.
       PARAMETER (NZBASELAT = -((59.3 D0/60 + 17)/60 +  41))    !Degrees South, thus negative.
       PARAMETER (NZBASELONG = ((34.65D0/60 + 46)/60 + 174))    !Degrees East.
C                                   Seconds  Minutes Degrees.
C   This is the location of the Mt. Cook trigonometrical base point for New Zealand.
C    (It's in the foyer of what used to be the Dominion Museum, Wellington)
 
       TYPE Terroir !Collate the attributes of location, as so far needed.
        CHARACTER*28 PLACENAME  !Name the location.
        DOUBLE PRECISION LATITUDE,LONGITUDE !Locate the location.
        INTEGER ZONETIME,ZONELONG   !The time zone of its civil clock.
       END TYPE Terroir     !The nature of the climate, soil, etc. is not yet involved.
       TYPE(Terroir) BASE   !Righto, let's have one of them.
       DATA BASE/Terroir("Mt. Cook Trig, Wellington.",  !The compiler bungles if NZBASEPLACE is used here.
     1  NZBASELAT,NZBASELONG,+12,180)/  !Where it's at.
 
      CONTAINS          !Let the madness begin.
 
       INTEGER*4 FUNCTION DAYNUM(YY,M,D)    !Computes JDayN - JDayShift, not JDayN.
C   Conversion from a Gregorian calendar date to a Julian day number, JDayN.
C   Valid for any Gregorian calendar date producing a Julian day number
C greater than zero, though remember that the Gregorian calendar
C was not used before y1582m10d15 and often, not after that either.
C   The zero of the Julian day number corresponds to the first of January
C 4713BC on the *Julian* calendar's naming scheme, as extended backwards
C into epochs when it did not exist. This function employs the naming
C scheme of the *Gregorian* calendar, and if extended backwards into
c epochs when it did not exist, it would compute a zero for y-4713m11d24
C not that this name was in use at the time...
C   Computationally, year 1 is preceeded by year 0, in a smooth progression.
C But there was never a year zero despite what astronomers like to say,
C so year 0 corresponds to 1BC, year -1 to 2BC, and so on back.
C Thus y-4713 in this counting would be 4714BC on the Gregorian calendar,
C were it to have existed then which it didn't.
C   To conform to the civil usage, the incoming YY, presumed a proper BC (negative)
C and AD (positive) year is converted into the computational counting sequence, Y,
C and used in the formula. If a YY = 0 is (improperly) offered, it will manifest
C as 1AD. Thus YY = -4714 will lead to calculations with Y = -4713.
C   For their convenience, astronomers decreed that a day starts at noon, so that
C in Europe, observations through the night all have the same day number.
C The current Western civil calendar however has the day starting just after midnight
C and that day's number lasts until the following midnight.
C
C   Example: Y = 1970, M = 1, D = 1;  JDAYN = 2440588, a Thursday but MOD(2440588,7) = 3.
C   and with the adjustment JDAYSHIFT, DAYNUM = 25568; mod 7 = 4 and DAYNAME(4) = "Thursday".
C
C   DAYNUM and MUNYAD are the infamous routine of H. F. Fliegel and T.C. van Flandern,
C   presented in Communications of the ACM, Vol. 11, No. 10 (October, 1968).
Carefully typed in by R.N.McLean (whom God preserve) December XXMMIIX.
C   Though I remain puzzled as to why they used I,J,K for Y,M,D,
C given that the variables were named in the INTEGER statement anyway.
        INTEGER*4 JDAYN     !Without rebasing, this won't fit in INTEGER*2.
        INTEGER YY,Y,M,D    !NB! Full year number, so 1970, not 70.
Caution: integer division in Fortran does not produce fractional results.
C The fractional part is discarded so that 4/3 gives 1 and -4/3 gives -1.
C Thus 4/3 might be Trunc(4/3) or 4 div 3 in other languages. Beware of negative numbers!
         Y = YY     !I can fiddle this copy without damaging the original's value.
         IF (Y.LT.1) Y = Y + 1  !Thus YY = -2=2BC, -1=1BC, +1=1AD becomes Y = -1, 0, 1..
         JDAYN = D - 32075
     a    + 1461*(Y + 4800  + (M - 14)/12)/4
     b    +  367*(M - 2     - (M - 14)/12*12)/12
     c    -    3*((Y + 4900 + (M - 14)/12)/100)/4
         DAYNUM = JDAYN - JDAYSHIFT !Thus, *NOT* the actual *Julian* Day Number.
       END FUNCTION DAYNUM      !But one such that Mod(n,7) gives day names.
 
Could compute the day of the year somewhat as follows...
c  DN:=D + (61*Month + (Month div 8)) div 2 - 30
c        + if Month > 2 then FebLength - 30 else 0;
 
       TYPE(DATEBAG) FUNCTION MUNYAD(DAYNUM)    !Oh for palindromic programming!
Conversion from a Julian day number to a Gregorian calendar date. See JDAYN/DAYNUM.
        INTEGER*4 DAYNUM,JDAYN  !Without rebasing, this won't fit in INTEGER*2.
        INTEGER Y,M,D,L,N       !Y will be a full year number: 1950 not 50.
         JDAYN = DAYNUM + JDAYSHIFT !Revert to a proper Julian day number.
         L = JDAYN + 68569  !Further machinations of H. F. Fliegel and T.C. van Flandern.
         N = 4*L/146097
         L = L - (146097*N + 3)/4
         Y = 4000*(L + 1)/1461001
         L = L - 1461*Y/4 + 31
         M = 80*L/2447
         D = L - 2447*M/80
         L = M/11
         M = M + 2 - 12*L
         Y = 100*(N - 49) + Y + L
         IF (Y.LT.1) Y = Y - 1  !The other side of conformity to BC/AD, as in DAYNUM.
         MUNYAD%YEAR  = Y   !Now place for the world to see.
         MUNYAD%MONTH = M
         MUNYAD%DAY   = D
       END FUNCTION MUNYAD  !A year has 365.2421988 days...
 
       INTEGER*4 FUNCTION EASTER(YEAR)  !Easter Sunday for the nominated year.
C    Apparently first given by an anonymous correspondent from New York to Nature in 1876.
C Samuel Butcher, Bishop of Meath, showed that this algorithm followed from Delambre's
C analytical solutions, and produces the date of Easter for all (Gregorian?) years.
C    The Julian and Gregorian calendars have abandoned the attempt to reconcile the
C period of the moon (a "moonth") with the length of the year, but Easter is a remnant
C of the time when moon-based calendars were preferred, and is a fertility rite aligned
C with the (northern hemisphere) spring, thus the Easter bunny and Easter eggs, and
C is named from the Mesopotamian goddess Ishtar.
C    Easter Sunday is the first Sunday after the first full moon on or after the vernal
c equinox of the northern hemisphere in March, which seems definite enough until details
c start to be considered. For instance, at what longitude is this to be decided
c from? Jerusalem? Like, the moment of full-moonness occurs at twenty-four
c different hours around the globe's time zones, not all of which are always
c in the same day. Similarly, the equinox may not be on the 21'st of March.
C    Whereupon, Easter Sunday is the Sunday following the Paschal Full Moon date for
c the year, and this is not quite what the actual moon does.
C    In June 325 A.D. astronomers approximated astronomical full moon dates for the
C Christian church, calling them Ecclesiastical Full Moon dates. From 326 A.D.
c the Paschal Full Moon date has always been the Ecclesiastical Full Moon date
c that followed the 20'th of March, which in 325AD was the equinox date.
C    The ecclesiastical Full Moon is defined as the fourteenth day of a tabular
C lunation, where day 1 corresponds to the ecclesiastical New Moon. The tables
c are based on the Metonic cycle, in which 235 mean synodic months occur in
c 6939.688 days. Since nineteen Gregorian years is 6939.6075 days, the dates
c of Moon phases in a given year will recur on nearly the same dates nineteen
c years later. To prevent the 0.08 day difference between the cycles from
c accumulating, the tables incorporate adjustments to synchronize the system
c over longer periods of time. Additional complications arise because the tabular
C lunations are of 29 or 30 integral days, the moon inconveniently failing to
c complete an orbit in a round number of days. The astronomical moon, as observed,
c might easily start a new moon before midnight in one location (local time)
c and after midnight in another, and so a nominal moon position is used,
c based on tables adjusted for "embolismic" lunations, and the head spins.
C    The entire system comprises a period of 5,700,000 years of 2,081,882,250 days,
c which is equated to 70,499,183 lunations. After this period, the dates of Easter
c repeat themselves.
C    Except, by then the length of the day and the period of the moon will have changed...
C
C    Variations on this routine are common. Be careful if the MOD operation might
C be presented with negative numbers, and on some computers the results will not
C be as expected. This routine does not involve any negative numbers, at least
C for positive years. Remember that this rule applies only from 326AD and that
C the Julian calendar was out of step with the seasons, leading to the introduction
C of the Gregorian calendar. As well, different churches defined Easter differently.
C
C    Earliest dates: 22/3/1818 (and 2285)
C    Latest dates:   25/4/1886 (and 1943, and 2038) for example.
 
       INTEGER a,Century,yy,d,e,f,g,h,i,k,l,m,p
       INTEGER Year,EasterMonth,EasterSunday
 
        a = MOD(year,19)        !a + 1 = Golden Number within the Metonic cycle.
        Century = year/100
        yy = MOD(year,100)
        d = Century/4;
        e = MOD(Century,4)
        f = (Century + 8)/25
        g = (Century - f + 1)/3
        h = MOD(19*a + Century - d - g + 15,30)
        i = yy/4
        k = MOD(yy,4)
        l = MOD(32 + 2*e + 2*i - h - k,7)   !32 + etc. prevents negative numbers.
        m = (a + 11*h + 22*l)/451
        EasterMonth = (h + l - 7*m + 114)/31;
        p = MOD(h + l - 7*m + 114,31)
        EasterSunday = p + 1            !Day in Easter Month.
        Easter = DAYNUM(Year,EasterMonth,EasterSunday)
       END FUNCTION EASTER
 
       DOUBLE PRECISION FUNCTION SOLARALT(EPOCH)    !Compute the Solar altitude, in degrees.
C   Location details are in the global variable BASE, with components LATITUDE, etc.
C   Entry is with a day number in terms of the civil calendar, as supplied by DAYNUM
c from a Gregorian Calendar date. The fractional part is the fraction of the day,
c so that noon corresponds to 0.5 of the day. To convert this to a Julian Day number
c with fractional part as used by astronomers requires adding back the JDAYSHIFT that
c was subtracted by DAYNUM, and stepping back half a day, since the civil day starts
c twelve hours before the astronomer's day.
c   And then remember that New Zealand is not on the longitude of Greenwich but twelve
c hours in advance, so that the time named by EPOCH in NZ at longitude ZONELONG is
c at Greenwich the moment named EPOCH - ZONETIME/24 in terms of days and fractions of days.
C Adjustments for daylight saving are for the caller to collate: the opinion
C here is that there is no daylight saving fiddle, so that noon is the solar noon, etc.
C   Confusions abound: angles in radians, degrees, minutes of a day, hours of a day, fraction of a day.
C   The calculations are based on page C24 of the Ephemeris for 2005, which offers low precision
c formulae for the Sun's coordinates and the equation of time. These are declared adequate
c to a hundreth of a degree of arc and a tenth of a minute in time between 1950 and 2050,
c and are centred on JDAY = 2451545.0, the Julian Day number and time as noon passed
c on 1/1/2000, at Greenwich. Remember that the solar disc is about half a degree across,
c and that atmospheric refraction raises its altitude by about half a degree at the horizon
c on average (but it can be double that): this calculation ignores such refraction.
C   Single precision (REAL*4) is not adequate for fine fractions of a day,
c since the day number is likely to be in the tens of thousands, thus consuming
c most of the six+ digits of precision available.
 
        DOUBLE PRECISION EPOCH  !Day and fraction of a day, local time (starting at midnight)...
        INTEGER*4 DAY       !A plain day number.
        DOUBLE PRECISION LOCALHOUR  !Within the day.
        DOUBLE PRECISION PI !A famous number.
        DATA PI/1D0/    !A "double precision" 1, rather than an approximate one.
        SAVE PI     !And keep it warm.
        DOUBLE PRECISION DEC,EOT,SOLARTIME,HA,T1,T2 !Time twiddling.
        DOUBLE PRECISION EN,EL,G,ELONG,EOBL,RA      !Time - Longitude exchanges, etc.
         IF (PI.EQ.1D0) PI = 4*ATAN(PI) !A notorious constant, evaluated in double precision.
         DAY = EPOCH            !Drop any fractional part of the day.
         LOCALHOUR = (EPOCH - DAY)*24   !The hour time of day (measured from midnight).
         EN = EPOCH - BASE%ZONETIME/24.0 + (JDAYSHIFT - 2451545) - 0.5  !Deviation from 1/1/2000, astronomically.
         EL = 280.460 + 0.9856474D0*EN  !Mean longitude of the sun, corrected for aberration..
         G =  357.528 + 0.9856003D0*EN  !Mean "anomaly".
         EL = EL - INT(EL/360)*360  !Adjust to [0:360).
         IF (EL.LT.0) EL = EL + 360 !The MOD function is unreliable for neg. values.
         G = G - INT(G/360)*360     !Adjust to [0:360)
         IF (G.LT.0) G = G + 360        !So, explicit working here...
         ELONG = EL + 1.915*SIND(G) + 0.020*SIND(2*G)   !Ecliptic longitude.
     ELONG = ELONG - INT(ELONG/360)*360 !These conversions are maddening.
         IF (ELONG.LT.0) ELONG = ELONG + 360    !+-180 rather than 0:360 might be better.
         EOBL = 23.439 - 0.0000004*EN       !Obliquity of ecliptic.
         RA = ATAN(COSD(EOBL)*TAND(ELONG))*180/PI   !Right ascension, in degrees.
         RA = RA - INT(RA/360)*360      !Expansion about a middle point,
         IF (RA.LT.0) RA = RA + 360     !Rather than one end of the range might be better too.
         IF (RA - ELONG.GT.180) RA = RA - 180   !Same quadrant as ELONG.
         IF (ELONG - RA.GT.180) RA = RA + 180   !So ensure this. I hope.
         DEC = ASIN(SIND(EOBL)*SIND(ELONG)) !Declination, in radians.
         EOT = (EL - RA)            !Equation of Time, in Degrees.
         IF (EOT.GT.+180) EOT = EOT - 360   !This is a difference between two directions.
         IF (EOT.LT.-180) EOT = EOT + 360   !And would be better as +-180.
         SOLARTIME = LOCALHOUR +        !Local civil time of day to...
     1    (EOT + BASE%LONGITUDE - BASE%ZONELONG)/15 !Sun's angle time, hours of a day.
         HA = (15*(SOLARTIME - 12))*PI/180  !Hour angle away from noon. Radians!
         T1 = SIN(DEC)*SIND(BASE%LATITUDE)  !Orientation of the pole towards the sun.
         T2 = COS(DEC)*COSD(BASE%LATITUDE)*COS(HA)  !As the Earth spins.
         SOLARALT = ASIN(T1 + T2)*180/PI    !Ah! Spherical trigonometry! (Apologies to Goethe)
       END FUNCTION SOLARALT    !Sun worship.
 
       CHARACTER*10 FUNCTION SLASHDATE(DAYNUM)  !This is relatively innocent.
Caution! The Gregorian calendar did not exist prior to 15/10/1582!
Confound fixed fields: because of the minus sign, years prior to 999BC will produce "****".
        INTEGER*4 DAYNUM    !-32768 to 32767 is just not adequate.
        TYPE(DATEBAG) D     !Though these numbers are more restrained.
         D = MUNYAD(DAYNUM) !Get the pieces.
         WRITE (SLASHDATE,1) D%DAY,D%MONTH,D%YEAR   !Some compilers will bungle this.
    1    FORMAT (I2,"/",I2,"/",I4)  !If so, a local variable must be used.
        RETURN              !As when SLASHDATE is invoked in a WRITE statement.
       END FUNCTION SLASHDATE   !Simple enough.
 
       CHARACTER*20 FUNCTION SAYDATE(D) !Unpack to civil form.
Caution! If D is negative, MOD may behave unhelpfully, thus the +7 and second MOD.
        INTEGER*4 D !The day.
         SAYDATE = DAYNAME(MOD(MOD(D,7) + 7,7))//" "//SLASHDATE(D)  !Take that.
       END FUNCTION SAYDATE !Ah, verbiage.
 
       CHARACTER*44 FUNCTION SAYDATES(IST,LST)  !One date good, two dates good good.
        INTEGER*4 IST,LST   !The day numbers.
        INTEGER L,LSTNB     !A finger.
        CHARACTER*44 SPLOT  !A work area.
         SPLOT = SAYDATE(IST)   !The start.
         L = LSTNB(SPLOT) + 1   !I'm not counting.
         IF (IST.NE.LST) SPLOT(L:) = " to " // SAYDATE(LST) !A proper span?
         SAYDATES = SPLOT   !Take that.
       END FUNCTION SAYDATES    !Ah, verbiage.
 
       CHARACTER*15 FUNCTION HMS(T) !Report the time of day.
Careful! Finite precision and binary/decimal/sexigesimal conversion could cause 2:30:00am. to appear as 2:29:60am.
        DOUBLE PRECISION S,T    !Seconds (completed) into the day.
        INTEGER H,M     !More traditional units are to be extracted.
        CHARACTER*15 TEXT   !A scratchpad.
         IF (T.EQ.SECONDSINDAY/2) THEN  !This might happen.
           TEXT = "High Noon!"  !Though the chances are thin.
         ELSE IF (T.EQ.SECONDSINDAY) THEN   !This is synonymous with the start of the next day.
           TEXT = "Midnight!"   !So this presumably won't happen.
         ELSE       !But more likely are miscellaneous values.
           H = T/3600       !Convert seconds into whole hours completed.
           S = T - H*3600   !The remaining time.
           M = S/60     !Seconds into minutes completed.
           S = S - M*60     !Remove them.
           IF (S .GT. 59.99949999999990D0) THEN !Likely to round up to 60?
             S = 0      !Yes. Curse recurring binary sequences for decimal.
             M = M + 1      !So, up the minute count.
             IF (M.GE.60) THEN  !Is there an overflow here too?
               M = 0        !Yes.
               H = H + 1    !So might appear 24:00:00.000 though not Midnight!
             END IF     !So much for twiddling the minutes.
           END IF       !And twiddling the hours.
           IF (H.LT.12) THEN    !A plague on the machine mentality.
             WRITE (TEXT,1) H,M,S,"am." !Ante-meridian.
    1        FORMAT (I2,":",I2,":",F6.3,A3) !Thus.
            ELSE        !For the post-meridian, H >= 12.
             IF (H.GT.12) H = H - 12    !Adjust to civil usage. NB! 12 appears.
             WRITE (TEXT,1) H,M,S,"pm." !Thus. Post-meridian.
           END IF   !So much for those fiddles.
           IF (TEXT(4:4).EQ." ") TEXT(4:4) = "0"    !Now help hint that the
           IF (TEXT(7:7).EQ." ") TEXT(7:7) = "0"    ! character string is one entity.
         END IF     !So much for preparation.
         HMS = TEXT !The result.
       END FUNCTION HMS !Possible compiler confusion if HMS is invoked in a WRITE statement.
 
       CHARACTER*39 FUNCTION TIMESTAMP(T)   !Unpack a time stamp.
        DOUBLE PRECISION T,C    !Full fractions.
        INTEGER*4 D !Integral day number.
         D = T      !Acquire the day number part, truncation NOT rounding.
         C = (T - D)*SECONDSINDAY   !The fraction of the day.
         TIMESTAMP = HMS(C)//" on "//SAYDATE(D) !Get it together.
       END FUNCTION TIMESTAMP   !So much for that.
 
       DOUBLE PRECISION FUNCTION NOWWAS(Idummy) !Ascertain the local time for interval assessment.
Compute with whole day numbers, to avoid day rollover annoyances.
Can't use single precision and expect much discrimination within a day.
C I'd prefer a TIMESTAMP(Local) and a TIMESTAMP(GMT) system function.
C Quite likely, the system separates its data to deliver the parameters, which I then re-glue.
        INTEGER IDummy  !A pacifier for the compiler.
        INTEGER MARK(8) !The computer's clock time will appear here, but fragmented.
         MARK(1) = IDummy   !Soothe the compiler by using the parameter.
         CALL DATE_AND_TIME(VALUES = MARK)  !Unpack info that I will repack.
c         WRITE (MSG,1) MARK
c    1    FORMAT ("The computer clock system reports:",
c     1    /"Year",I5,", Month",I3,", Day",I3,
c     2    /" Minutes from GMT",I5,
c     3    /" Hour",I3,", Minute",I3,",Seconds",I3,".",I3)
         NOWWAS = DAYNUM(MARK(1),MARK(2),MARK(3))   !The days stream past...
     1   + ((MARK(5)*60 + MARK(6))*60 + MARK(7) + MARK(8)/1000D0)   !By the millisecond.
     2     /SECONDSINDAY    !Fraction of a day, always less than 1 as MARK(5) is declared < 24.
       END FUNCTION NOWWAS  !The Hand of Time has already moved on.
 
       SUBROUTINE NZDAYJOLT(YEAR,STRETCH,SHRINK)    !Daylight saving changeover days...
Compute the transition days, in this YEAR, even though a summer time span crosses a year boundary.
Can't be bothered with the original half-hour shift of NZ's time zone introduced during WW1.
        INTEGER YEAR,STRETCH,SHRINK !From the year, determine the trouble.
        INTEGER Y,I         !Assistants.
Calculate the start of summer time for this YEAR.
         Y = MAX(YEAR,1975)     !The first daylight saving start was in 1975.
         IF (Y < 1989) THEN     !This was the early rule.
           I = DAYNUM(Y,10,31)      !Find the end of October.
           I = I - MOD(I,7)     !Thus the last Sunday in October.
         ELSE IF (Y.EQ.1989) THEN   !A transition.
           I = DAYNUM(1989,10,8)    !A specified date. (The second Sunday)
         ELSE               !And now the new rule.
           I = DAYNUM(Y,10,1)       !Find the first of October.
           I = I + MOD(7 - MOD(I,7),7)  !Thus the first Sunday.
         END IF             !Possibly, yet further variations?
         SHRINK = I         !Pass the word.
Calculate the end of summer time for YEAR. STRETCH *precedes* SHRINK in the southern hemisphere.
         Y = MAX(YEAR,1976) !The first summertime ending was in 1976.
         I = DAYNUM(Y,3,1)      !The first of March.
         I = I + MOD(7 - MOD(I,7),7)    !Advance to Sunday, if not one.
         IF (Y > 1989) I = I + 14   !Thus the third Sunday.
         STRETCH = I            !This day has 25 hours in it.
       END SUBROUTINE NZDAYJOLT !Others shall confront the collateral damage.

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c Subroutine by: D. Bala Subrahamanyam (subrahamanyam@gmail.com)       !
c Last modified: (ATJ/DBS: March 24, 2010)                             !
c This subroutine will count the day number starting from the 1st day  !
c dates (e.g., between 13.11.1979 to 22.01.2010) ...                   !
c Input to subroutine:                                                 !
c               ----> First Date     (idd)                             !
c               ----> First Month    (imm)                             !
c               ----> First Year     (iyy)                             !
c               ----> Second Date    (jdd)                             !         
c               ----> Second Month   (jmm)                             !
c               ----> Second Year    (jyr)                             !        
c Output of the surbroutine:                                           ! 
c               ----> No. of days between first & second date (nday)  !
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
        Subroutine makedates (idd,imm,iyy,jdd,jmm,jyy,nod)
 
        Call caljul (idd, imm, iyy, jday1)
        Call caljul (jdd, jmm, jyy, jday2)
 
        If (jyy .lt. iyy) Write (*, *) 'WRONG INPUTS...'
        If (jyy .lt. iyy) Write (*, *) 'FIRST DATE > SECOND DATE'
        If (jyy .lt. iyy) Go to 200
        
 
        Call caljul (31, 12, iyy, lday1)
        Call caljul (31, 12, jyy, lday2)
 
        If ((jyy .eq. iyy).and.(lday1 .gt. lday2)) Write (*, *) 'WRONG INPUTS' 
        If ((jyy .eq. iyy).and.(lday1 .gt lday2)) Write (*, *) 'WRONG INPUTS' 
 
 
        nod1 = lday1 - jday1
        nod2 = jday2
 
        Do 10 i = iyy+1, jyy-1
 
        Call caljul (31, 12, i, jadd)
 
        iadd = iadd + jadd
 
 10     continue
 
        If (jyy .eq. iyy) nod = (jday2 - jday1) + 1
        If (jyy .eq. iyy+1) nod = nod1 + nod2 + 1
        If (jyy .gt. iyy+1) nod = nod1 + iadd + nod2 + 1                 
 
 200    Return
 
        End
  
 
 
        Subroutine caljul (iddd, immm, iyyy, julday)
 
        Dimension incmon(13), leapinc(13)
 
        Data incmon /0, 31, 59, 90, 120, 151, 181, 212,   
     &               243, 273, 304, 334, 365/
 
        Data leapinc /0, 31, 60, 91, 121, 152, 182, 213,  
     &               244, 274, 305, 335, 366/
 
        icheck = mod (iyyy, 4) 
c       icheck is leap year checking index
c       icheck = 0, if the current year is not a leap year
c       icheck = 1, if the current year is a leap year
 
        if (icheck .ne. 0) julday = iddd + incmon(immm)
        if (icheck .eq. 0) julday = iddd + leapinc(immm)
       
        Return
 
        End