! 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
!####################################################################
!
!
! 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).
!
!
! 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).
!
!
! call orbit
!
!
!
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
!###################################################################
!
!
! 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.
!
!
! 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.
!
!
! r = r_inv_squared (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
!
!
!
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
!####################################################################
!
!
! 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.
!
!
! 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.
!
!
! r = angle (t)
!
!
! time of year (between 0 and 2*pi; t=0 at NH autumnal
! equinox
!
!
!
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
!####################################################################
!
!
! declination returns the solar declination angle at orbital
! position ang in earth's orbit.
!
!
! declination returns the solar declination angle at orbital
! position ang in earth's orbit.
!
!
! r = declination (ang)
!
!
! solar orbital position ang in earth's orbit
!
!
!
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