! A program to compute the sun declination for each day of the year.
!
! Functionality copied from astronomy.F90

program dailydeclination

  implicit none
  
! use

  real, parameter :: PI      = 3.14159265358979323846
  real, parameter :: RAD_TO_DEG=180./PI
  real, parameter :: DEG_TO_RAD=PI/180.

  real, parameter :: seconds_per_day=86400.   ! seconds in a day
  real, parameter :: twopi = 2*PI

  !-------- namelist  ---------

  real    :: ecc   = 0.01671   ! eccentricity of orbital ellipse 
                               ! [ non-dimensional ]
  real    :: obliq = 23.439    ! obliquity [ degrees ]
  real    :: per   = 102.932   ! longitude of perihelion with respect 
                               ! to autumnal equinox in NH [ degrees ]
  !integer :: period = 0        ! specified length of year [ seconds ] ; 
  !                             ! must be specified to override default 
  !                             ! value given by length_of_year in 
  !                             ! time_manager_mod
  !integer :: day_ae    = 23    ! day of specified autumnal equinox
  !integer :: month_ae  = 9     ! month of specified autumnal equinox
  !integer :: year_ae   = 1998  ! year of specified autumnal equinox
  !integer :: hour_ae   = 5     ! hour of specified autumnal equinox
  !integer :: minute_ae = 37    ! minute of specified autumnal equinox
  !integer :: second_ae = 0     ! second of specified autumnal equinox
  integer :: num_angles = 3600 ! number of intervals into which the year 
                               ! is divided to compute orbital positions

  real, dimension(:), allocatable ::      &
       orb_angle ! table of orbital positions (0 to 
                 ! 2*pi) as a function of time  used
                 ! to find actual orbital position
                 ! via interpolation

  !   time of year; autumnal equinox = 0.0,
  !                    one year = 2 * pi
  real :: time_since_ae, ang, dec

  integer :: day_of_year, equinox_day, day_since_equinox


  equinox_day = 31+28+31+30+31+30+31+31+23 ! 23rd september

  !---------------------------------------------------------------------
  !    call orbit to define table of orbital angles as function of 
  !    orbital time.
  !----------------------------------------------------------------------
  ! wfc moved here from orbit
  allocate ( orb_angle(0:num_angles) )
  call orbit


  !---------------------------------------------------------------------
  !    define the orbital angle (location in year), solar declination and
  !    earth sun distance factor. use functions contained in this module.
  !---------------------------------------------------------------------
  do day_of_year=1,365 ! day of year
     day_since_equinox = abs(mod(day_of_year-1 + (365-equinox_day+1), 365))
     !write(*,*) day_of_year, equinox_day, day_since_equinox
     time_since_ae = day_since_equinox/365.0*360.0*deg_to_rad
     !--------------------------------------------------------------------
     !    be sure the time in the annual cycle is legitimate.
     !---------------------------------------------------------------------
     if (time_since_ae < 0.0 .or. time_since_ae > twopi) then
        write(*,*) 'time_since_ae not between 0 and 2pi', time_since_ae
        stop
     end if
     ang = angle(time_since_ae)
     dec = declination(ang)
     write(*,'(i3,a,f8.5,a,f8.5,a,f8.5,a,f6.2)') &
          day_of_year,' ',time_since_ae,' ',ang,' ',dec,' ',dec*rad_to_deg
  end do

contains

!####################################################################
! <SUBROUTINE NAME="orbit">
!  <OVERVIEW>
!    orbit computes and stores a table of value of orbital angles as a 
!    function of orbital time (both the angle and time are zero at 
!    autumnal equinox in the NH, and range from 0 to 2*pi).
!  </OVERVIEW>
!  <DESCRIPTION>
!    orbit computes and stores a table of value of orbital angles as a 
!    function of orbital time (both the angle and time are zero at 
!    autumnal equinox in the NH, and range from 0 to 2*pi).
!  </DESCRIPTION>
!  <TEMPLATE>
!   call orbit
!  </TEMPLATE>
! </SUBROUTINE>
!
subroutine orbit

!---------------------------------------------------------------------
!    orbit computes and stores a table of value of orbital angles as a 
!    function of orbital time (both the angle and time are zero at 
!    autumnal equinox in the NH, and range from 0 to 2*pi).
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!   local variables

      integer :: n
      real    :: d1, d2, d3, d4, d5, dt, norm

!--------------------------------------------------------------------
!    allocate the orbital angle array, sized by the namelist parameter
!    num_angles, defining the annual cycle resolution of the earth's
!    orbit. define some constants to be used.
!--------------------------------------------------------------------
! wfc moving to astronomy_init
!     allocate ( orb_angle(0:num_angles) )
      orb_angle(0) = 0.0
      dt = twopi/float(num_angles)
      norm = sqrt(1.0 - ecc**2)
      dt = dt*norm

!---------------------------------------------------------------------
!    define the orbital angle at each of the num_angles locations in 
!    the orbit.
!---------------------------------------------------------------------
      do n = 1, num_angles
        d1 = dt*r_inv_squared(orb_angle(n-1))
        d2 = dt*r_inv_squared(orb_angle(n-1)+0.5*d1)
        d3 = dt*r_inv_squared(orb_angle(n-1)+0.5*d2)
        d4 = dt*r_inv_squared(orb_angle(n-1)+d3)
        d5 = d1/6.0 + d2/3.0 +d3/3.0 +d4/6.0
        orb_angle(n) = orb_angle(n-1) + d5
      end do
  
!-------------------------------------------------------------------



end subroutine orbit


!###################################################################
! <FUNCTION NAME="r_inv_squared">
!  <OVERVIEW>
!    r_inv_squared returns the inverse of the square of the earth-sun
!    distance relative to the mean distance at angle ang in the earth's
!    orbit.
!  </OVERVIEW>
!  <DESCRIPTION>
!    r_inv_squared returns the inverse of the square of the earth-sun
!    distance relative to the mean distance at angle ang in the earth's
!    orbit.
!  </DESCRIPTION>
!  <TEMPLATE>
!    r = r_inv_squared (ang)
!  </TEMPLATE>
!  <IN NAME="ang" TYPE="real">
!    angular position of earth in its orbit, relative to a 
!            value of 0.0 at the NH autumnal equinox, value between
!            0.0 and 2 * pi
!  </IN>
! </FUNCTION>
!
function r_inv_squared (ang)

!--------------------------------------------------------------------
!    r_inv_squared returns the inverse of the square of the earth-sun
!    distance relative to the mean distance at angle ang in the earth's
!    orbit.
!--------------------------------------------------------------------

!--------------------------------------------------------------------
real, intent(in) :: ang
!--------------------------------------------------------------------

!---------------------------------------------------------------------
!
!  intent(in) variables:
!
!      ang   angular position of earth in its orbit, relative to a 
!            value of 0.0 at the NH autumnal equinox, value between
!            0.0 and 2 * pi
!            [ radians ]
!
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!  local variables
 
      real :: r_inv_squared, r, rad_per

!---------------------------------------------------------------------
!  local variables:
!
!      r_inv_squared    the inverse of the square of the earth-sun
!                       distance relative to the mean distance 
!                       [ dimensionless ]
!      r                earth-sun distance relative to mean distance
!                       [ dimensionless ]
!      rad_per          angular position of perihelion 
!                       [ radians ]
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!    define the earth-sun distance (r) and then return the inverse of
!    its square (r_inv_squared) to the calling routine.
!--------------------------------------------------------------------
      rad_per       = per*deg_to_rad
      r             = (1 - ecc**2)/(1. + ecc*cos(ang - rad_per))
      r_inv_squared = r**(-2)


end function r_inv_squared


!####################################################################
! <FUNCTION NAME="angle">
!  <OVERVIEW>
!    angle determines the position within the earth's orbit at time t
!    in the year ( t = 0 at NH autumnal equinox.) by interpolating
!    into the orbital position table. 
!  </OVERVIEW>
!  <DESCRIPTION>
!    angle determines the position within the earth's orbit at time t
!    in the year ( t = 0 at NH autumnal equinox.) by interpolating
!    into the orbital position table. 
!  </DESCRIPTION>
!  <TEMPLATE>
!    r = angle (t)
!  </TEMPLATE>
!  <IN NAME="t" TYPE="real">
!    time of year (between 0 and 2*pi; t=0 at NH autumnal
!                equinox
!  </IN>
! </FUNCTION>
!
function angle (t)

!--------------------------------------------------------------------
!    angle determines the position within the earth's orbit at time t
!    in the year ( t = 0 at NH autumnal equinox.) by interpolating
!    into the orbital position table. 
!--------------------------------------------------------------------

!--------------------------------------------------------------------
real, intent(in) :: t
!--------------------------------------------------------------------

!-------------------------------------------------------------------
!
!  intent(in) variables:
!
!         t      time of year (between 0 and 2*pi; t=0 at NH autumnal
!                equinox
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables

      real :: angle, norm_time, x
      integer :: int, int_1

!--------------------------------------------------------------------
!   local variables:
!
!     angle       orbital position relative to NH autumnal equinox
!                 [ radians ]
!     norm_time   index into orbital table corresponding to input time
!                 [ dimensionless ]
!     x           fractional distance between the orbital table entries
!                 bracketing the input time
!                 [ dimensionless ]
!     int         table index which is lower than actual position, but
!                 closest to it
!                 [ dimensionless ]
!     int_1       next table index just larger than actual orbital 
!                 position
!                 [ dimensionless ]
!
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!    define orbital tables indices bracketing current orbital time
!    (int and int_1). define table index distance between the lower 
!    table value (int) and the actual orbital time (x). define orbital
!    position as being  x of the way between int and int_1. renormalize
!    angle to be within the range 0 to 2*pi.
!--------------------------------------------------------------------
      norm_time = t*float(num_angles)/twopi
      int = floor(norm_time)
      int = modulo(int,num_angles)
      int_1 = int+1
      x = norm_time - floor(norm_time)
      angle = (1.0 -x)*orb_angle(int) + x*orb_angle(int_1)
      angle = modulo(angle, twopi)

end function angle


!####################################################################
! <FUNCTION NAME="declination">
!  <OVERVIEW>
!    declination returns the solar declination angle at orbital
!    position ang in earth's orbit.
!  </OVERVIEW>
!  <DESCRIPTION>
!    declination returns the solar declination angle at orbital
!    position ang in earth's orbit.
!  </DESCRIPTION>
!  <TEMPLATE>
!    r =  declination (ang)
!  </TEMPLATE>
!  <IN NAME="ang" TYPE="real">
!     solar orbital position ang in earth's orbit
!  </IN>
! </FUNCTION>
!
function declination (ang)

!--------------------------------------------------------------------
!    declination returns the solar declination angle at orbital
!    position ang in earth's orbit.
!--------------------------------------------------------------------

!--------------------------------------------------------------------
real, intent(in) :: ang
!--------------------------------------------------------------------

!--------------------------------------------------------------------
!   local variables

      real :: declination
      real :: rad_obliq, sin_dec

!--------------------------------------------------------------------
!   local variables:
!
!    declination         solar declination angle
!                        [ radians ]              
!    rad_obliq           obliquity of the ecliptic
!                        [ radians ]              
!    sin_dec             sine of the solar declination
!                        [ dimensionless ]
!
!--------------------------------------------------------------------

!---------------------------------------------------------------------
!    compute the solar declination.
!---------------------------------------------------------------------
      rad_obliq   =   obliq*deg_to_rad
      sin_dec     = - sin(rad_obliq)*sin(ang)
      declination =   asin(sin_dec)


end function declination

end program dailydeclination
