I am looking for feedback on a couple of functions/subroutines for converting between Julian and Gregorian calendar in Fortran. It is for an internal software for extreme value analysis, so it needs to be robust.
I am currently using something I found online (credit: http://aa.usno.navy.mil/faq/docs/JD_Formula.php) with a slight modification to suit my needs of also having a fractional time added to the Julian date.
To convert a date vector (Gregorian format) to Julian and fractional time:
INTEGER FUNCTION JULIAN(YEAR,MONTH,DAY)
IMPLICIT NONE
INTEGER, INTENT(IN) :: YEAR,MONTH,DAY
JULIAN = DAY-32075+1461*(YEAR+4800+(MONTH-14)/12)/4+367*(MONTH-2-(MONTH-14)/12*12)/12-3*((YEAR+4900+(MONTH-14)/12)/100)/4
END FUNCTION
REAL(DOUBLE) FUNCTION FRACTIONTIME(HOUR,MINUTE,SECOND)
IMPLICIT NONE
INTEGER, INTENT(IN) :: HOUR,MINUTE,SECOND
FRACTIONTIME = (HOUR + (MINUTE + SECOND/60.D0)/60.D0)/24.D0
END FUNCTION FRACTIONTIME
The Julian date I store is then JD = Julian(...) + Fractiontime(...)
And for converting back to Gregorian format:
SUBROUTINE GREGORIAN(JD,YEAR,MONTH,DAY,HOUR,MINUTE,SECOND)
IMPLICIT NONE
REAL(DOUBLE), INTENT(IN) :: JD
INTEGER, INTENT(OUT) :: YEAR, MONTH, DAY, HOUR, MINUTE, SECOND
REAL(DOUBLE) :: JT
INTEGER :: I,J,K,L,N
L = INT(JD)+68569
N = 4*L/146097
L = L-(146097*N+3)/4
I = 4000*(L+1)/1461001
L = L-1461*I/4+31
J = 80*L/2447
K = L-2447*J/80
L = J/11
J = J+2-12*L
I = 100*(N-49)+I+L
YEAR = I
MONTH = J
DAY = K
JT = DMOD(JD,1.D0)*24.D0
HOUR = INT(JT)
JT = DMOD(JT,1.D0)*60.D0
MINUTE = INT(JT)
JT = DMOD(JT,1.D0)*60.D0
SECOND = NINT(JT)
IF (SECOND == 60) THEN
SECOND = SECOND-60
MINUTE = MINUTE+1
END IF
END SUBROUTINE GREGORIAN
If this is perfect as-is, then see this as a way to share the work, but if you see any problems or have suggestions for other ways to tackle this, I am all ears!
Short description of what I am actually doing:
- Time series of e.g. wind speeds are stored in a standard format with a date vector (year, month, day, hour, minute, second) as the first 6 columns and then the wind speed in another.
- Peaks in the wind speeds are identified, which need to be statistically independent, i.e. more than 72 hours apart
- In order to ensure this, I store the Julian and then test for any neighbouring peaks that are JD+3 and JD-3 apart
- For output reasons and because I don't want to drag my original date vector around, I convert the Julian back to Gregorian
And to keep in line with the FAQ:
What is the best way to convert between Julian and Gregorian in Fortran 95/2003?
