! bof
! **********************************************************************
! Fortran 95 module to provide julian date computation

! **********************************************************************
! Source Control Strings

! $Id: julian.f90 1.3 2001/06/01 23:07:26Z Dan Release $

! **********************************************************************
! Copyright 2000 Purple Sage Computing Soltuions, Inc.

! **********************************************************************
! Julian dates allow the number of days between dates to be computed

! **********************************************************************
! Summary of License

!   This library is free software; you can redistribute it and/or
!   modify it under the terms of the GNU Library General Public
!   License as published by the Free Software Foundation; either
!   version 2 of the License, or (at your option) any later version.

!   This library is distributed in the hope that it will be useful,
!   but WITHOUT ANY WARRANTY; without even the implied warranty of
!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!   Library General Public License for more details.

!   You should have received a copy of the GNU Library General Public
!   License along with this library; if not, write to the Free
!   Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

! To report bugs, suggest enhancements, etc. to the Authors,
! Contact:
!    Purple Sage Computing Solutions, Inc.
!                               send email to dnagle@@erols.com
!                                   or fax to 703 471 0684 (USA)
!                                  or mail to 12142 Purple Sage Ct.
!                                             Reston, VA 20194-5621 USA

! *********************************************************************

!  Module to provide the calculation of Julian dates,
!  after the C routines found in _ACM SIGNUM Newsletter_
!  v 30, n 4, October 1995, beginning on p 14,
!  "Algorithms for Julian Dates" by Richard L. Branham, Jr.

!  Dates B.C. should be expressed in astronomical format,
!  that is, date B.C. = -(date-1).  Thus, 747 B.C. = -746.
!  Days are expressd as per the astronomical convention
!  of days beginning and ending at noon rather than the
!  civil convention of midnight.  The astronomical day
!  is equal to the civil day minus 0.5.  Thus, the civil
!  15th of the month is the astronomical day 15.5.

!  Translated to Fortran 90 by Dan Nagle, January, 1996.

!  Modifications by Dan Nagle
!  Use array of cumulative days per month from beginning
!     of year to replace the loops with nested conditional
!     code found in the original C code
!  julian_day_of_week() has day as double precision argument in
!     conformity with the other routines rather than int
!     as in the original C code

!  Routines added by Dan Nagle, Febuary, 2000.
!  julian_now()
!  julian_today()
!  julian_fract()
!  julian_hms()
!  julian_astro_to_civil()
!  julian_civil_to_astro()
!  julian_solar()
!  julian_solar_ymd()

! *********************************************************************

!  julian_date library

!     julian() julian date of a month, day, year

!     julian_now() julian date of now
!     julian_today() julian date of noon today

!     julian_ymd() year, month, day of a julian date

!     julian_day_of_week() character string given month, day, year

!     julian_fract() fraction of a day given hour, minute, second
!     julian_hms() hour, minute, second given julian fraction

!     julian_astro_to_civil() change convention
!     julian_civil_to_astro()

!     julian_solar() eighth of year, day of eighth, given year, month, day
!     julian_solar_ymd() year, month, day given year,  eighth of year, day of eighth

! *********************************************************************

!  julian_date

module julian_date

! *********************************************************************

!  julian_date uses

!  import double_k

use standard_types

!  import .mod.

use standard_functions

! *********************************************************************

!  force explicit declarations

implicit none                                                        ! declare all

! *********************************************************************

!  force explicit exports

private                                                              ! export all

! *********************************************************************

!  force declarations

character( len= *), public, parameter :: julian_rcs_id = &
   '$Id: julian.f90 1.3 2001/06/01 23:07:26Z Dan Release $'

! *********************************************************************

!  days prior to beginning of each month (non&) leap year

integer, parameter, dimension( 13) :: &
   bgn_day = (/ 0,  31,  59,  90, 120, 151, &
              181, 212, 243, 273, 304, 334, 365 /), &
   bgn_day_ly = (/ 0,  31,  60,  91, 121, 152, &
                 182, 213, 244, 274, 305, 335, 366 /)

! *********************************************************************

!  julian_date library

public :: julian

public :: julian_now
public :: julian_today

public :: julian_ymd

public :: julian_day_of_week

public :: julian_fract
public :: julian_hms

public :: julian_astro_to_civil
public :: julian_civil_to_astro

public :: julian_solar
public :: julian_solar_ymd

! *********************************************************************

!  module procedures

contains                                                             ! julian_date

!  Function calculates a Julian date for a given year,
!  month, day on the Julian calendar.  Dates B.C.
!  should be expressed in astronomical format, that is,
!  date B.C. = -(date-1).  Thus, 747 B.C. = -746.
!  Days are expressd as per the astronomical convention
!  of days beginning and ending at noon rather than the
!  civil convention of midnight.  The astronomical day
!  is equal to the civil day minus 0.5

!  Input to the function julian are the year and month
!  expressed as integers, day as a double precision;
!  output is the Julian date expressed as a double precision

! *********************************************************************

!  julian()

pure real( kind= double_k) function julian( year, month, day)

! *********************************************************************

!  julian() interface

integer, intent( in) :: year                                         ! starting 4713 BC

integer, intent( in) :: month                                        ! 1 thru 12

real( kind= double_k), intent( in) :: day                            ! civil day

! *********************************************************************

!  julian() local

   real( kind= double_k) :: astro_day                                ! starts 12 hr before civil day

! *********************************************************************

!  julian() text

continue                                                             ! julian()

!  correct for difference between astronomical and civil time

   astro_day = day - 0.5_double_k                                    ! civil day less 12 hours

!  calculate julian day for beginning of year (1461 days in 4 years starting 4713 BC)

   julian = real( (( 1461*( 4712 + year) - 1) / 4), kind= double_k)

!  1 Jan 4713 BC is Julian Date Zero

   origin_year: if( year == -4712 )then

      julian = julian - 1.0_double_k                                 ! -4712 is a leap year

   endif origin_year

!  add in number of days in current month

   julian = julian + astro_day                                       ! year + day

!  add in number of days in preceeding months leap year (&non)

   leap_year: if( ( year .mod. 4) == 0 )then                         ! leap year or not

      julian = julian + real( bgn_day_ly( month), double_k)

   else leap_year                                                    ! leap year or not

      julian = julian + real( bgn_day( month), double_k)

   endif leap_year                                                   ! leap year or not

!  julian date

return                                                               ! julian()

! *********************************************************************

!  julian()

end function julian

! *********************************************************************

!  Subroutine to take a Julian date and return a year,
!  month, day on the Julian calendar.  The caller
!  provides the integers year, month and the double precision
!  day.  julian_ymd returns double precision jd.

! *********************************************************************

!  julian_ymd()

pure subroutine julian_ymd( jd, year, month, day)

! *********************************************************************

!  julian_ymd() interface

real( kind= double_k), intent( in) ::  jd                            ! julian date

integer, intent( out) :: year                                        ! X BC == -(X-1)

integer, intent( out) :: month                                       ! 1 ==> 12

real( kind= double_k), intent( out) ::  day                          ! civil day

! *********************************************************************

!  julian_ymd() text

continue                                                             ! julian_ymd()

!  julian date zero is a special case

   jd_near_zero: if( abs( jd) <= epsilon( jd) )then                  ! jd == 0.0

      year = -4712                                                   ! -4712 ==> 4713 BC
      month = 1                                                      ! 1 ==> January
      day = 0.5_double_k                                             ! civil day 1 ==> jd 0.5

      return                                                         ! julian_ymd()

   endif jd_near_zero                                                ! jd == 0.0

!  number of julian years in julian date less origin of julian system

   year = int( ( jd / 365.25_double_k) - 4712.0_double_k)            ! integer year

!  correct for X BC == -(X-1)

   bc_year: if( jd <= 0.0_double_k )then

      year = year - 1                                                ! BC convention

   endif bc_year

!  number of days since beginning of year

   day = jd - real( ( 1461*( 4712 + year) - 1) / 4, kind= double_k) + 0.5_double_k

!  calculate month

   month = 1                                                         ! start with January

!  (julian) leap years special case

   month_leap_year: if( ( year .mod. 4) == 0 )then                   ! leap year or not

      ly_month: do while ( int( day) > bgn_day_ly( month + 1) )

         month = month + 1                                           ! while day > days in month so far

      enddo ly_month

      day = day - real( bgn_day_ly( month), kind= double_k)

   else month_leap_year                                              ! leap year or not

      nly_month: do while ( int( day) > bgn_day( month + 1) )

         month = month + 1                                           ! while day > days in month so far

      enddo nly_month

      day = day - real( bgn_day( month), kind= double_k)

   endif month_leap_year                                             ! leap year or not

!  julian_ymd()

return                                                               ! julian_ymd()

! *********************************************************************

!  julian_ymd()

end subroutine julian_ymd

! *********************************************************************

!  Subroutine to calculate the day of the week corresponding
!  to the year, month, day on the Julian calendar.  The
!  function returns its results in the character variable dow.

! *********************************************************************

!  julian_day_of_week()

pure subroutine julian_day_of_week( dow, year, month, day)

! *********************************************************************

!  julian_day_of_week() interface

integer, intent( in) :: year                                         ! year of date

integer, intent( in) :: month                                        ! month of date

real( kind= double_k), intent( in) :: day                            ! day of date

character( len= *), intent( out) :: dow                              ! day of week

! *********************************************************************

!  julian_day_of_week() local

   integer :: yr                                                     ! local year
   integer :: mon                                                    ! local month
   integer :: id                                                     ! local day

   integer :: n                                                      ! jd modulo seven

! *********************************************************************

!  julian_day_of_week() text

continue                                                             ! julian_day_of_week()

   id = nint( day)                                                   ! nearest whole day

!  all years positive (by adding multiple of 28)

   yr = year + 4732                                                  ! 4732 > 4713

!  treat January and February as months 13 & 14 of preceeding year

   jan_feb: if( month <= 2 )then                                     ! adjust January & Febuary

      mon = month + 12                                               ! months 1 & 2 ==> 13 & 14
      yr = yr - 1                                                    ! of the preceeding year

   else jan_feb                                                      ! adjust January & Febuary

      mon = month                                                    ! month is month

   endif jan_feb                                                     ! adjust January & Febuary

!  calculate numerical parameter (herein, magic numbers: see Branham)

   n = ( ( id + 2*mon + ( 3*mon + 3 )/5 + yr + yr/4 + 6) .mod. 7) + 1

!  select dow name string or null

   get_string: select case( n)                                       ! 1 <= n <= 7

   case( 1) get_string

      dow = 'Monday'                                                 ! Luna

   case( 2) get_string

      dow = 'Tuesday'                                                ! Mars

   case( 3) get_string

      dow = 'Wednesday'                                              ! Mercury

   case( 4) get_string

      dow = 'Thursday'                                               ! Jupiter

   case( 5) get_string

      dow = 'Friday'                                                 ! Venus

   case( 6) get_string

      dow = 'Saturday'                                               ! Saturn

   case( 7) get_string

      dow = 'Sunday'                                                 ! Sol

   case default get_string

      dow = ''                                                       ! null string

   end select get_string                                             ! 1 <= n <= 7

!  dow set

return                                                               ! julian_day_of_week()

! *********************************************************************

!  julian_day_of_week()

end subroutine julian_day_of_week

! *********************************************************************

!  julian_today()

pure real( kind= double_k) function julian_today()

! *********************************************************************

!  julian_today() local

   integer, dimension( values_size) :: v_array

! *********************************************************************

!  julian_today() text

continue                                                             ! julian_today()

!  get year, month, day of today

   call date_and_time( values= v_array)

!  calculate julian day for beginning of year ( using magic numbers)

   julian_today = real( (( 1461*( 4712 + v_array( values_year)) - 1) / 4), double_k)

!  add in number of days in current month

   julian_today = julian_today + real( v_array( values_day), kind= double_k)

!  civil to astro

   julian_today = julian_today - 0.5_double_k

!  add in number of days in preceeding months leap year (&non)

   leap_year: if( ( v_array( values_year) .mod. 4) == 0 )then

      julian_today = julian_today + real( bgn_day_ly( v_array( values_month)), double_k)

   else leap_year

      julian_today = julian_today + real( bgn_day( v_array( values_month)), double_k)

   endif leap_year

!  julian_today()

return                                                               ! julian_today()

! *********************************************************************

!  julian_today()

end function julian_today

! *********************************************************************

!  julian_now()

pure real( kind= double_k) function julian_now()

! *********************************************************************

!  julian_now() local

   integer, dimension( values_size) :: v_array

   integer :: cnow, cmax

! *********************************************************************

!  julian_now() text

continue                                                             ! julian_now()

!  get year, month, day of today

   call date_and_time( values= v_array)

!  calculate julian day for beginning of year ( using magic numbers)

   julian_now = real( (( 1461*( 4712 + v_array( values_year)) - 1) / 4), double_k)

!  add in number of days in current month

   julian_now = julian_now + real( v_array( values_day), kind= double_k)

!  add fraction of today

   call system_clock( count= cnow, count_max= cmax)

   julian_now = julian_now + real( cnow, kind= double_k) / real( cmax, kind= double_k)

!  civil to astro

   julian_now = julian_now - 0.5_double_k

!  add in number of days in preceeding months leap year (&non)

   leap_year: if( ( v_array( values_year) .mod. 4) == 0 )then

      julian_now = julian_now + real( bgn_day_ly( v_array( values_month)), double_k)

   else leap_year

      julian_now = julian_now + real( bgn_day( v_array( values_month)), double_k)

   endif leap_year

!  julian_now()

return                                                               ! julian_now()

! *********************************************************************

!  julian_now()

end function julian_now

! *********************************************************************

!  julian_fract()

real( kind= double_k) function julian_fract( hour, minute, second, msec)

! *********************************************************************

!  julian_fract() interface

integer, intent( in) :: hour, minute, second, msec

! *********************************************************************

!  julian_fract() assumes a 24 hour clock

! *********************************************************************

!  julian_fract() text

continue                                                             ! julian_fract()

!  compute real seconds

   julian_fract = real( second, kind= double_k) + real( msec, kind= double_k)/1000.0_double_k

!  compute real minutes

   julian_fract = real( minute, kind= double_k) + julian_fract/60.0_double_k

!  compute real hours

   julian_fract = real( hour, kind= double_k) + julian_fract/60.0_double_k

!  julian fraction of a day

   julian_fract = julian_fract/24.0_double_k

return                                                               ! julian_fract()

! *********************************************************************

!  julian_fract()

end function julian_fract

! *********************************************************************

!  julian_hms()

pure subroutine julian_hms( jd, hour, minute, second, msec)

! *********************************************************************

!  julian_hms() interface

real( kind= double_k), intent( in) :: jd

integer, intent( out) :: hour, minute, second, msec

! *********************************************************************

!  julian_hms() local

   real( kind= double_k) :: fract

! *********************************************************************

!  julian_hms() assumes a 24 hour clock

! *********************************************************************

!  julian_hms() text

continue                                                             ! julian_hms()

!  compute hour

   hour = int( jd * 24.0_double_k)                                   ! 24 hours/day

   fract = ( jd - real( hour, kind= double_k) / 24.0_double_k) * 24.0_double_k

!  compute minute

   minute = int( fract * 60.0_double_k)                              ! 60 minutes/hour

   fract = ( fract - real( minute, kind= double_k) / 60.0_double_k) * 60.0_double_k

!  compute second

   second = int( fract * 60.0_double_k)                              ! 60 seconds/minute

   fract = ( fract - real( second, kind= double_k) / 60.0_double_k) * 60.0_double_k

!  compute msec

   msec = nint( fract * 1000.0_double_k)                             ! 1000 msec/second

return                                                               ! julian_hms()

! *********************************************************************

!  julian_hms()

end subroutine julian_hms

! *********************************************************************

!  julian_astro_to_civil()

pure real( kind= double_k) function julian_astro_to_civil( astro_day)

! *********************************************************************

!  julian_astro_to_civil() interface

real( kind= double_k), intent( in) :: astro_day

! *********************************************************************

!  julian_astro_to_civil() text

continue                                                             ! julian_astro_to_civil()

!  correct for difference between astronomical and civil time

   julian_astro_to_civil = astro_day + 0.5_double_k

return                                                               ! julian_astro_to_civil()

! *********************************************************************

!  julian_astro_to_civil()

end function julian_astro_to_civil

! *********************************************************************

!  julian_civil_to_astro()

pure real( kind= double_k) function julian_civil_to_astro( civil_day)

! *********************************************************************

!  julian_civil_to_astro() interface

real( kind= double_k), intent( in) :: civil_day

! *********************************************************************

!  julian_civil_to_astro() text

continue                                                             ! julian_civil_to_astro()

!  correct for difference between astronomical and civil time

   julian_civil_to_astro = civil_day - 0.5_double_k

return                                                               ! julian_civil_to_astro()

! *********************************************************************

!  julian_civil_to_astro()

end function julian_civil_to_astro

! *********************************************************************

!  julian_solar()

pure subroutine julian_solar( eoy, doe, year, month, day)

! *********************************************************************

!  julian_solar() interface

integer, intent( out) :: eoy

integer, intent( out) :: doe

integer, intent( inout) :: year

integer, intent( in) :: month

real( kind= double_k), intent( in) ::  day

! *********************************************************************

!  julian_solar() local

   integer :: iday, idelta

! *********************************************************************

!  julian_solar() text

continue                                                             ! julian_solar()

   iday = int( day + 0.5_double_k)                                   ! astro to civil

   if( month == 12 .and. iday >= 21 )then                            ! after winter solstice

      idelta = int( julian( year, month, day) - julian( year, month, real( 21, kind= double_k)) ) + 1

      year = year - 1

   else                                                              ! before winter solstice

      idelta = int( julian( year, month, day) - julian( year - 1, 12, real( 21, kind= double_k)) ) + 1

   endif                                                             ! before or after winter solstice

   eoy = idelta / 8

   doe = idelta - eoy*8

return                                                               ! julian_solar()

! *********************************************************************

!  julian_solar()

end subroutine julian_solar

! *********************************************************************

!  julian_solar_ymd()

pure subroutine julian_solar_ymd( eoy, doe, year, month, day)

! *********************************************************************

!  julian_solar_ymd() interface

integer, intent( in) :: eoy

integer, intent( in) :: doe

integer, intent( inout) :: year

integer, intent( out) :: month

real( kind= double_k), intent( out) :: day

! *********************************************************************

!  julian_solar_ymd() local

   real( kind= double_k), parameter :: dpe = 365.25_double_k / 8.0_double_k

   integer :: ieoy

   real( kind= double_k) :: jd

! *********************************************************************

!  julian_solar_ymd() text

continue                                                             ! julian_solar_ymd()

   if( eoy == 1 .and. doe <= 11 )then

      jd = julian( year - 1, 12, 20.5_double_k + real( doe, kind= double_k))

      call julian_ymd( jd, year, month, day)

   else

      ieoy = eoy - 1

      jd = julian( year - 1, 12, 20.5_double_k) + real( nint( ieoy*dpe) + doe, kind= double_k)

      call julian_ymd( jd, year, month, day)

   endif

return                                                               ! julian_solar_ymd()

! *********************************************************************

!  julian_solar_ymd()

end subroutine julian_solar_ymd

! **********************************************************************

!  julian_date

! **********************************************************************
! $Id: julian.f90 1.3 2001/06/01 23:07:26Z Dan Release $

end module julian_date                                               ! eof