!===============================================================================
! CVS: $Id: shr_orb_mod.F90,v 1.2.4.1 2004/01/02 18:50:55 mvr Exp $
! CVS: $Source: /fs/cgd/csm/models/CVS.REPOS/shared/csm_share/shr/shr_orb_mod.F90,v $
! CVS: $Name: cam3_0_brnchT_release01 $
!===============================================================================
MODULE shr_orb_mod 4,3
use shr_kind_mod
, only: SHR_KIND_R8, SHR_KIND_IN
use shr_sys_mod, only: shr_sys_abort
use shr_const_mod, only: SHR_CONST_PI
IMPLICIT none
!----------------------------------------------------------------------------
! PUBLIC: Interfaces and global data
!----------------------------------------------------------------------------
public shr_orb_cosz, shr_orb_params, shr_orb_decl, shr_orb_print
real (SHR_KIND_R8),public, parameter :: SHR_ORB_UNDEF_REAL = 1.e36 ! undefined real
integer(SHR_KIND_IN),public, parameter :: SHR_ORB_UNDEF_INT = 2000000000 ! undefined int
!----------------------------------------------------------------------------
! PRIVATE: by default everything else is private to this module
!----------------------------------------------------------------------------
private
real (SHR_KIND_R8),parameter :: pi = SHR_CONST_PI
real (SHR_KIND_R8),parameter :: SHR_ORB_ECCEN_MIN = 0.0 ! min value for eccen
real (SHR_KIND_R8),parameter :: SHR_ORB_ECCEN_MAX = 0.1 ! max value for eccen
real (SHR_KIND_R8),parameter :: SHR_ORB_OBLIQ_MIN = -90.0 ! min value for obliq
real (SHR_KIND_R8),parameter :: SHR_ORB_OBLIQ_MAX = +90.0 ! max value for obliq
real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MIN = 0.0 ! min value for mvelp
real (SHR_KIND_R8),parameter :: SHR_ORB_MVELP_MAX = 360.0 ! max value for mvelp
CONTAINS
!===============================================================================
real(SHR_KIND_R8) FUNCTION shr_orb_cosz(jday,lat,lon,declin) 1
!----------------------------------------------------------------------------
!
! FUNCTION to return the cosine of the solar zenith angle.
! Assumes 365.0 days/year.
!
!--------------- Code History -----------------------------------------------
!
! Original Author: Brian Kauffman
! Date: Jan/98
! History: adapted from statement FUNCTION in share/orb_cosz.h
!
!----------------------------------------------------------------------------
real (SHR_KIND_R8),intent(in) :: jday ! Julian cal day (1.xx to 365.xx)
real (SHR_KIND_R8),intent(in) :: lat ! Centered latitude (radians)
real (SHR_KIND_R8),intent(in) :: lon ! Centered longitude (radians)
real (SHR_KIND_R8),intent(in) :: declin ! Solar declination (radians)
!----------------------------------------------------------------------------
shr_orb_cosz = sin(lat)*sin(declin) - &
& cos(lat)*cos(declin)*cos(jday*2.0*pi + lon)
END FUNCTION shr_orb_cosz
!===============================================================================
SUBROUTINE shr_orb_params( iyear_AD , eccen , obliq , mvelp , & 1,5
& obliqr , lambm0 , mvelpp, log_print )
!-------------------------------------------------------------------------------
!
! Calculate earths orbital parameters using Dave Threshers formula which
! came from Berger, Andre. 1978 "A Simple Algorithm to Compute Long-Term
! Variations of Daily Insolation". Contribution 18, Institute of Astronomy
! and Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium
!
!------------------------------Code history-------------------------------------
!
! Original Author: Erik Kluzek
! Date: Oct/97
!
!-------------------------------------------------------------------------------
!----------------------------- Arguments ------------------------------------
real (SHR_KIND_R8),intent(inout) :: eccen ! orbital eccentricity
real (SHR_KIND_R8),intent(inout) :: obliq ! obliquity in degrees
real (SHR_KIND_R8),intent(inout) :: mvelp ! moving vernal equinox long
integer(SHR_KIND_IN),intent(in) :: iyear_AD ! Year to calculate orbit for
logical ,intent(in) :: log_print ! Flags print of status/error
real (SHR_KIND_R8),intent(out) :: obliqr ! Earths obliquity in rad
real (SHR_KIND_R8),intent(out) :: lambm0 ! Mean long of perihelion at
! vernal equinox (radians)
real (SHR_KIND_R8),intent(out) :: mvelpp ! moving vernal equinox long
! of perihelion plus pi (rad)
!------------------------------ Parameters ----------------------------------
integer(SHR_KIND_IN),parameter :: poblen =47 ! # of elements in series wrt obliquity
integer(SHR_KIND_IN),parameter :: pecclen=19 ! # of elements in series wrt eccentricity
integer(SHR_KIND_IN),parameter :: pmvelen=78 ! # of elements in series wrt vernal equinox
real (SHR_KIND_R8),parameter :: psecdeg = 1.0/3600.0 ! arc sec to deg conversion
real (SHR_KIND_R8) :: degrad = pi/180. ! degree to radian conversion factor
real (SHR_KIND_R8) :: yb4_1950AD ! number of years before 1950 AD
! Cosine series data for computation of obliquity: amplitude (arc seconds),
! rate (arc seconds/year), phase (degrees).
real (SHR_KIND_R8), parameter :: obamp(poblen) = & ! amplitudes for obliquity cos series
& (/ -2462.2214466, -857.3232075, -629.3231835, &
& -414.2804924, -311.7632587, 308.9408604, &
& -162.5533601, -116.1077911, 101.1189923, &
& -67.6856209, 24.9079067, 22.5811241, &
& -21.1648355, -15.6549876, 15.3936813, &
& 14.6660938, -11.7273029, 10.2742696, &
& 6.4914588, 5.8539148, -5.4872205, &
& -5.4290191, 5.1609570, 5.0786314, &
& -4.0735782, 3.7227167, 3.3971932, &
& -2.8347004, -2.6550721, -2.5717867, &
& -2.4712188, 2.4625410, 2.2464112, &
& -2.0755511, -1.9713669, -1.8813061, &
& -1.8468785, 1.8186742, 1.7601888, &
& -1.5428851, 1.4738838, -1.4593669, &
& 1.4192259, -1.1818980, 1.1756474, &
& -1.1316126, 1.0896928/)
real (SHR_KIND_R8), parameter :: obrate(poblen) = & ! rates for obliquity cosine series
& (/ 31.609974, 32.620504, 24.172203, &
& 31.983787, 44.828336, 30.973257, &
& 43.668246, 32.246691, 30.599444, &
& 42.681324, 43.836462, 47.439436, &
& 63.219948, 64.230478, 1.010530, &
& 7.437771, 55.782177, 0.373813, &
& 13.218362, 62.583231, 63.593761, &
& 76.438310, 45.815258, 8.448301, &
& 56.792707, 49.747842, 12.058272, &
& 75.278220, 65.241008, 64.604291, &
& 1.647247, 7.811584, 12.207832, &
& 63.856665, 56.155990, 77.448840, &
& 6.801054, 62.209418, 20.656133, &
& 48.344406, 55.145460, 69.000539, &
& 11.071350, 74.291298, 11.047742, &
& 0.636717, 12.844549/)
real (SHR_KIND_R8), parameter :: obphas(poblen) = & ! phases for obliquity cosine series
& (/ 251.9025, 280.8325, 128.3057, &
& 292.7252, 15.3747, 263.7951, &
& 308.4258, 240.0099, 222.9725, &
& 268.7809, 316.7998, 319.6024, &
& 143.8050, 172.7351, 28.9300, &
& 123.5968, 20.2082, 40.8226, &
& 123.4722, 155.6977, 184.6277, &
& 267.2772, 55.0196, 152.5268, &
& 49.1382, 204.6609, 56.5233, &
& 200.3284, 201.6651, 213.5577, &
& 17.0374, 164.4194, 94.5422, &
& 131.9124, 61.0309, 296.2073, &
& 135.4894, 114.8750, 247.0691, &
& 256.6114, 32.1008, 143.6804, &
& 16.8784, 160.6835, 27.5932, &
& 348.1074, 82.6496/)
! Cosine/sine series data for computation of eccentricity and fixed vernal
! equinox longitude of perihelion (fvelp): amplitude,
! rate (arc seconds/year), phase (degrees).
real (SHR_KIND_R8), parameter :: ecamp (pecclen) = & ! ampl for eccen/fvelp cos/sin series
& (/ 0.01860798, 0.01627522, -0.01300660, &
& 0.00988829, -0.00336700, 0.00333077, &
& -0.00235400, 0.00140015, 0.00100700, &
& 0.00085700, 0.00064990, 0.00059900, &
& 0.00037800, -0.00033700, 0.00027600, &
& 0.00018200, -0.00017400, -0.00012400, &
& 0.00001250/)
real (SHR_KIND_R8), parameter :: ecrate(pecclen) = & ! rates for eccen/fvelp cos/sin series
& (/ 4.2072050, 7.3460910, 17.8572630, &
& 17.2205460, 16.8467330, 5.1990790, &
& 18.2310760, 26.2167580, 6.3591690, &
& 16.2100160, 3.0651810, 16.5838290, &
& 18.4939800, 6.1909530, 18.8677930, &
& 17.4255670, 6.1860010, 18.4174410, &
& 0.6678630/)
real (SHR_KIND_R8), parameter :: ecphas(pecclen) = & ! phases for eccen/fvelp cos/sin series
& (/ 28.620089, 193.788772, 308.307024, &
& 320.199637, 279.376984, 87.195000, &
& 349.129677, 128.443387, 154.143880, &
& 291.269597, 114.860583, 332.092251, &
& 296.414411, 145.769910, 337.237063, &
& 152.092288, 126.839891, 210.667199, &
& 72.108838/)
! Sine series data for computation of moving vernal equinox longitude of
! perihelion: amplitude (arc seconds), rate (arc sec/year), phase (degrees).
real (SHR_KIND_R8), parameter :: mvamp (pmvelen) = & ! amplitudes for mvelp sine series
& (/ 7391.0225890, 2555.1526947, 2022.7629188, &
& -1973.6517951, 1240.2321818, 953.8679112, &
& -931.7537108, 872.3795383, 606.3544732, &
& -496.0274038, 456.9608039, 346.9462320, &
& -305.8412902, 249.6173246, -199.1027200, &
& 191.0560889, -175.2936572, 165.9068833, &
& 161.1285917, 139.7878093, -133.5228399, &
& 117.0673811, 104.6907281, 95.3227476, &
& 86.7824524, 86.0857729, 70.5893698, &
& -69.9719343, -62.5817473, 61.5450059, &
& -57.9364011, 57.1899832, -57.0236109, &
& -54.2119253, 53.2834147, 52.1223575, &
& -49.0059908, -48.3118757, -45.4191685, &
& -42.2357920, -34.7971099, 34.4623613, &
& -33.8356643, 33.6689362, -31.2521586, &
& -30.8798701, 28.4640769, -27.1960802, &
& 27.0860736, -26.3437456, 24.7253740, &
& 24.6732126, 24.4272733, 24.0127327, &
& 21.7150294, -21.5375347, 18.1148363, &
& -16.9603104, -16.1765215, 15.5567653, &
& 15.4846529, 15.2150632, 14.5047426, &
& -14.3873316, 13.1351419, 12.8776311, &
& 11.9867234, 11.9385578, 11.7030822, &
& 11.6018181, -11.2617293, -10.4664199, &
& 10.4333970, -10.2377466, 10.1934446, &
& -10.1280191, 10.0289441, -10.0034259/)
real (SHR_KIND_R8), parameter :: mvrate(pmvelen) = & ! rates for mvelp sine series
& (/ 31.609974, 32.620504, 24.172203, &
& 0.636717, 31.983787, 3.138886, &
& 30.973257, 44.828336, 0.991874, &
& 0.373813, 43.668246, 32.246691, &
& 30.599444, 2.147012, 10.511172, &
& 42.681324, 13.650058, 0.986922, &
& 9.874455, 13.013341, 0.262904, &
& 0.004952, 1.142024, 63.219948, &
& 0.205021, 2.151964, 64.230478, &
& 43.836462, 47.439436, 1.384343, &
& 7.437771, 18.829299, 9.500642, &
& 0.431696, 1.160090, 55.782177, &
& 12.639528, 1.155138, 0.168216, &
& 1.647247, 10.884985, 5.610937, &
& 12.658184, 1.010530, 1.983748, &
& 14.023871, 0.560178, 1.273434, &
& 12.021467, 62.583231, 63.593761, &
& 76.438310, 4.280910, 13.218362, &
& 17.818769, 8.359495, 56.792707, &
& 8.448301, 1.978796, 8.863925, &
& 0.186365, 8.996212, 6.771027, &
& 45.815258, 12.002811, 75.278220, &
& 65.241008, 18.870667, 22.009553, &
& 64.604291, 11.498094, 0.578834, &
& 9.237738, 49.747842, 2.147012, &
& 1.196895, 2.133898, 0.173168/)
real (SHR_KIND_R8), parameter :: mvphas(pmvelen) = & ! phases for mvelp sine series
& (/ 251.9025, 280.8325, 128.3057, &
& 348.1074, 292.7252, 165.1686, &
& 263.7951, 15.3747, 58.5749, &
& 40.8226, 308.4258, 240.0099, &
& 222.9725, 106.5937, 114.5182, &
& 268.7809, 279.6869, 39.6448, &
& 126.4108, 291.5795, 307.2848, &
& 18.9300, 273.7596, 143.8050, &
& 191.8927, 125.5237, 172.7351, &
& 316.7998, 319.6024, 69.7526, &
& 123.5968, 217.6432, 85.5882, &
& 156.2147, 66.9489, 20.2082, &
& 250.7568, 48.0188, 8.3739, &
& 17.0374, 155.3409, 94.1709, &
& 221.1120, 28.9300, 117.1498, &
& 320.5095, 262.3602, 336.2148, &
& 233.0046, 155.6977, 184.6277, &
& 267.2772, 78.9281, 123.4722, &
& 188.7132, 180.1364, 49.1382, &
& 152.5268, 98.2198, 97.4808, &
& 221.5376, 168.2438, 161.1199, &
& 55.0196, 262.6495, 200.3284, &
& 201.6651, 294.6547, 99.8233, &
& 213.5577, 154.1631, 232.7153, &
& 138.3034, 204.6609, 106.5938, &
& 250.4676, 332.3345, 27.3039/)
!---------------------------Local variables----------------------------------
integer(SHR_KIND_IN) :: i ! Index for series summations
real (SHR_KIND_R8) :: obsum ! Obliquity series summation
real (SHR_KIND_R8) :: cossum ! Cos series summation for eccentricity/fvelp
real (SHR_KIND_R8) :: sinsum ! Sin series summation for eccentricity/fvelp
real (SHR_KIND_R8) :: fvelp ! Fixed vernal equinox long of perihelion
real (SHR_KIND_R8) :: mvsum ! mvelp series summation
real (SHR_KIND_R8) :: beta ! Intermediate argument for lambm0
real (SHR_KIND_R8) :: years ! Years to time of interest ( pos <=> future)
real (SHR_KIND_R8) :: eccen2 ! eccentricity squared
real (SHR_KIND_R8) :: eccen3 ! eccentricity cubed
!-------------------------- Formats -----------------------------------------
character(len=*),parameter :: F00 = "('(shr_orb_params) ',4a)"
character(len=*),parameter :: F01 = "('(shr_orb_params) ',a,i9)"
character(len=*),parameter :: F02 = "('(shr_orb_params) ',a,f6.3)"
character(len=*),parameter :: F03 = "('(shr_orb_params) ',a,es14.6)"
!----------------------------------------------------------------------------
! radinp and algorithms below will need a degree to radian conversion factor
if ( log_print ) then
write(6,F00) 'Calculate characteristics of the orbit:'
write(6,F00) 'CVS revision: $Revision: 1.2.4.1 $'
write(6,F00) 'CVS Tag : $Name: cam3_0_brnchT_release01 $'
end if
! Check for flag to use input orbit parameters
IF ( iyear_AD == SHR_ORB_UNDEF_INT ) THEN
! Check input obliq, eccen, and mvelp to ensure reasonable
if( obliq == SHR_ORB_UNDEF_REAL )then
if ( log_print ) then
write(6,F00) 'Have to specify orbital parameters:'
write(6,F00) 'Either set: iyear_AD, OR [obliq, eccen, and mvelp]:'
write(6,F00) 'iyear_AD is the year to simulate orbit for (ie. 1950): '
write(6,F00) 'obliq, eccen, mvelp specify the orbit directly:'
write(6,F00) 'The AMIP II settings (for a 1995 orbit) are: '
write(6,F00) ' obliq = 23.4441'
write(6,F00) ' eccen = 0.016715'
write(6,F00) ' mvelp = 102.7'
end if
call shr_sys_abort()
else if ( log_print ) then
write(6,F00) 'Use input orbital parameters: '
end if
if( (obliq < SHR_ORB_OBLIQ_MIN).or.(obliq > SHR_ORB_OBLIQ_MAX) ) then
if ( log_print ) then
write(6,F03) 'Input obliquity unreasonable: ', obliq
end if
call shr_sys_abort()
end if
if( (eccen < SHR_ORB_ECCEN_MIN).or.(eccen > SHR_ORB_ECCEN_MAX) ) then
if ( log_print ) then
write(6,F03) 'Input eccentricity unreasonable: ', eccen
end if
call shr_sys_abort()
end if
if( (mvelp < SHR_ORB_MVELP_MIN).or.(mvelp > SHR_ORB_MVELP_MAX) ) then
if ( log_print ) then
write(6,F03) 'Input mvelp unreasonable: ' , mvelp
end if
call shr_sys_abort()
end if
eccen2 = eccen*eccen
eccen3 = eccen2*eccen
ELSE ! Otherwise calculate based on years before present
if ( log_print ) then
write(6,F01) 'Calculate orbit for year: ' , iyear_AD
end if
yb4_1950AD = 1950.0 - float(iyear_AD)
if ( abs(yb4_1950AD) .gt. 1000000.0 )then
if ( log_print ) then
write(6,F00) 'orbit only valid for years+-1000000'
write(6,F00) 'Relative to 1950 AD'
write(6,F03) '# of years before 1950: ',yb4_1950AD
write(6,F01) 'Year to simulate was : ',iyear_AD
end if
call shr_sys_abort()
end if
! The following calculates the earths obliquity, orbital eccentricity
! (and various powers of it) and vernal equinox mean longitude of
! perihelion for years in the past (future = negative of years past),
! using constants (see parameter section) given in the program of:
!
! Berger, Andre. 1978 A Simple Algorithm to Compute Long-Term Variations
! of Daily Insolation. Contribution 18, Institute of Astronomy and
! Geophysics, Universite Catholique de Louvain, Louvain-la-Neuve, Belgium.
!
! and formulas given in the paper (where less precise constants are also
! given):
!
! Berger, Andre. 1978. Long-Term Variations of Daily Insolation and
! Quaternary Climatic Changes. J. of the Atmo. Sci. 35:2362-2367
!
! The algorithm is valid only to 1,000,000 years past or hence.
! For a solution valid to 5-10 million years past see the above author.
! Algorithm below is better for years closer to present than is the
! 5-10 million year solution.
!
! Years to time of interest must be negative of years before present
! (1950) in formulas that follow.
years = - yb4_1950AD
! In the summations below, cosine or sine arguments, which end up in
! degrees, must be converted to radians via multiplication by degrad.
!
! Summation of cosine series for obliquity (epsilon in Berger 1978) in
! degrees. Convert the amplitudes and rates, which are in arc secs, into
! degrees via multiplication by psecdeg (arc seconds to degrees conversion
! factor). For obliq, first term is Berger 1978 epsilon star; second
! term is series summation in degrees.
obsum = 0.0
do i = 1, poblen
obsum = obsum + obamp(i)*psecdeg*cos((obrate(i)*psecdeg*years + &
& obphas(i))*degrad)
end do
obliq = 23.320556 + obsum
! Summation of cosine and sine series for computation of eccentricity
! (eccen; e in Berger 1978) and fixed vernal equinox longitude of
! perihelion (fvelp; pi in Berger 1978), which is used for computation
! of moving vernal equinox longitude of perihelion. Convert the rates,
! which are in arc seconds, into degrees via multiplication by psecdeg.
cossum = 0.0
do i = 1, pecclen
cossum = cossum+ecamp(i)*cos((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
end do
sinsum = 0.0
do i = 1, pecclen
sinsum = sinsum+ecamp(i)*sin((ecrate(i)*psecdeg*years+ecphas(i))*degrad)
end do
! Use summations to calculate eccentricity
eccen2 = cossum*cossum + sinsum*sinsum
eccen = sqrt(eccen2)
eccen3 = eccen2*eccen
! A series of cases for fvelp, which is in radians.
if (abs(cossum) .le. 1.0E-8) then
if (sinsum .eq. 0.0) then
fvelp = 0.0
else if (sinsum .lt. 0.0) then
fvelp = 1.5*pi
else if (sinsum .gt. 0.0) then
fvelp = .5*pi
endif
else if (cossum .lt. 0.0) then
fvelp = atan(sinsum/cossum) + pi
else if (cossum .gt. 0.0) then
if (sinsum .lt. 0.0) then
fvelp = atan(sinsum/cossum) + 2.0*pi
else
fvelp = atan(sinsum/cossum)
endif
endif
! Summation of sin series for computation of moving vernal equinox long
! of perihelion (mvelp; omega bar in Berger 1978) in degrees. For mvelp,
! first term is fvelp in degrees; second term is Berger 1978 psi bar
! times years and in degrees; third term is Berger 1978 zeta; fourth
! term is series summation in degrees. Convert the amplitudes and rates,
! which are in arc seconds, into degrees via multiplication by psecdeg.
! Series summation plus second and third terms constitute Berger 1978
! psi, which is the general precession.
mvsum = 0.0
do i = 1, pmvelen
mvsum = mvsum + mvamp(i)*psecdeg*sin((mvrate(i)*psecdeg*years + &
& mvphas(i))*degrad)
end do
mvelp = fvelp/degrad + 50.439273*psecdeg*years + 3.392506 + mvsum
! Cases to make sure mvelp is between 0 and 360.
do while (mvelp .lt. 0.0)
mvelp = mvelp + 360.0
end do
do while (mvelp .ge. 360.0)
mvelp = mvelp - 360.0
end do
END IF ! end of test on whether to calculate or use input orbital params
! Orbit needs the obliquity in radians
obliqr = obliq*degrad
! 180 degrees must be added to mvelp since observations are made from the
! earth and the sun is considered (wrongly for the algorithm) to go around
! the earth. For a more graphic explanation see Appendix B in:
!
! A. Berger, M. Loutre and C. Tricot. 1993. Insolation and Earth Orbital
! Periods. J. of Geophysical Research 98:10,341-10,362.
!
! Additionally, orbit will need this value in radians. So mvelp becomes
! mvelpp (mvelp plus pi)
mvelpp = (mvelp + 180.)*degrad
! Set up an argument used several times in lambm0 calculation ahead.
beta = sqrt(1. - eccen2)
! The mean longitude at the vernal equinox (lambda m nought in Berger
! 1978; in radians) is calculated from the following formula given in
! Berger 1978. At the vernal equinox the true longitude (lambda in Berger
! 1978) is 0.
lambm0 = 2.*((.5*eccen + .125*eccen3)*(1. + beta)*sin(mvelpp) &
& - .250*eccen2*(.5 + beta)*sin(2.*mvelpp) &
& + .125*eccen3*(1./3. + beta)*sin(3.*mvelpp))
if ( log_print ) then
write(6,F03) '------ Computed Orbital Parameters ------'
write(6,F03) 'Eccentricity = ',eccen
write(6,F03) 'Obliquity (deg) = ',obliq
write(6,F03) 'Obliquity (rad) = ',obliqr
write(6,F03) 'Long of perh(deg) = ',mvelp
write(6,F03) 'Long of perh(rad) = ',mvelpp
write(6,F03) 'Long at v.e.(rad) = ',lambm0
write(6,F03) '-----------------------------------------'
end if
END SUBROUTINE shr_orb_params
!===============================================================================
SUBROUTINE shr_orb_decl(calday ,eccen ,mvelpp ,lambm0 ,obliqr ,delta ,eccf) 2
!-------------------------------------------------------------------------------
!
! Compute earth/orbit parameters using formula suggested by
! Duane Thresher.
!
!---------------------------Code history----------------------------------------
!
! Original version: Erik Kluzek
! Date: Oct/1997
!
!-------------------------------------------------------------------------------
!------------------------------Arguments--------------------------------
real (SHR_KIND_R8),intent(in) :: calday ! Calendar day, including fraction
real (SHR_KIND_R8),intent(in) :: eccen ! Eccentricity
real (SHR_KIND_R8),intent(in) :: obliqr ! Earths obliquity in radians
real (SHR_KIND_R8),intent(in) :: lambm0 ! Mean long of perihelion at the
! vernal equinox (radians)
real (SHR_KIND_R8),intent(in) :: mvelpp ! moving vernal equinox longitude
! of perihelion plus pi (radians)
real (SHR_KIND_R8),intent(out) :: delta ! Solar declination angle in rad
real (SHR_KIND_R8),intent(out) :: eccf ! Earth-sun distance factor (ie. (1/r)**2)
!---------------------------Local variables-----------------------------
real (SHR_KIND_R8),parameter :: dayspy = 365.0 ! days per year
real (SHR_KIND_R8),parameter :: ve = 80.5 ! Calday of vernal equinox
! assumes Jan 1 = calday 1
real (SHR_KIND_R8) :: lambm ! Lambda m, mean long of perihelion (rad)
real (SHR_KIND_R8) :: lmm ! Intermediate argument involving lambm
real (SHR_KIND_R8) :: lamb ! Lambda, the earths long of perihelion
real (SHR_KIND_R8) :: invrho ! Inverse normalized sun/earth distance
real (SHR_KIND_R8) :: sinl ! Sine of lmm
! Compute eccentricity factor and solar declination using
! day value where a round day (such as 213.0) refers to 0z at
! Greenwich longitude.
!
! Use formulas from Berger, Andre 1978: Long-Term Variations of Daily
! Insolation and Quaternary Climatic Changes. J. of the Atmo. Sci.
! 35:2362-2367.
!
! To get the earths true longitude (position in orbit; lambda in Berger
! 1978) which is necessary to find the eccentricity factor and declination,
! must first calculate the mean longitude (lambda m in Berger 1978) at
! the present day. This is done by adding to lambm0 (the mean longitude
! at the vernal equinox, set as March 21 at noon, when lambda=0; in radians)
! an increment (delta lambda m in Berger 1978) that is the number of
! days past or before (a negative increment) the vernal equinox divided by
! the days in a model year times the 2*pi radians in a complete orbit.
lambm = lambm0 + (calday - ve)*2.*pi/dayspy
lmm = lambm - mvelpp
! The earths true longitude, in radians, is then found from
! the formula in Berger 1978:
sinl = sin(lmm)
lamb = lambm + eccen*(2.*sinl + eccen*(1.25*sin(2.*lmm) &
& + eccen*((13.0/12.0)*sin(3.*lmm) - 0.25*sinl)))
! Using the obliquity, eccentricity, moving vernal equinox longitude of
! perihelion (plus), and earths true longitude, the declination (delta)
! and the normalized earth/sun distance (rho in Berger 1978; actually inverse
! rho will be used), and thus the eccentricity factor (eccf), can be
! calculated from formulas given in Berger 1978.
invrho = (1. + eccen*cos(lamb - mvelpp)) / (1. - eccen*eccen)
! Set solar declination and eccentricity factor
delta = asin(sin(obliqr)*sin(lamb))
eccf = invrho*invrho
return
END SUBROUTINE shr_orb_decl
!===============================================================================
SUBROUTINE shr_orb_print( iyear_AD, eccen, obliq, mvelp ) 1
!-------------------------------------------------------------------------------
!
! Print out the information on the Earths input orbital characteristics
!
!---------------------------Code history----------------------------------------
!
! Original version: Erik Kluzek
! Date: Oct/1997
!
!-------------------------------------------------------------------------------
!---------------------------Arguments----------------------------------------
integer(SHR_KIND_IN),intent(in) :: iyear_AD ! requested Year (AD)
real (SHR_KIND_R8),intent(in) :: eccen ! eccentricity (unitless)
! (typically 0 to 0.1)
real (SHR_KIND_R8),intent(in) :: obliq ! obliquity (-90 to +90 degrees)
! typically 22-26
real (SHR_KIND_R8),intent(in) :: mvelp ! moving vernal equinox at perhel
! (0 to 360 degrees)
!-------------------------- Formats -----------------------------------------
character(len=*),parameter :: F00 = "('(shr_orb_print) ',4a)"
character(len=*),parameter :: F01 = "('(shr_orb_print) ',a,i9.4)"
character(len=*),parameter :: F02 = "('(shr_orb_print) ',a,f6.3)"
character(len=*),parameter :: F03 = "('(shr_orb_print) ',a,es14.6)"
!----------------------------------------------------------------------------
if ( iyear_AD .ne. SHR_ORB_UNDEF_INT ) then
if ( iyear_AD > 0 ) then
write(6,F01) 'Orbital parameters calculated for year: AD ',iyear_AD
else
write(6,F01) 'Orbital parameters calculated for year: BC ',iyear_AD
end if
else if ( obliq /= SHR_ORB_UNDEF_REAL ) then
write(6,F03) 'Orbital parameters: '
write(6,F03) 'Obliquity (degree): ', obliq
write(6,F03) 'Eccentricity (unitless): ', eccen
write(6,F03) 'Long. of moving Perhelion (deg): ', mvelp
else
write(6,F03) 'Orbit parameters not set!'
end if
END SUBROUTINE shr_orb_print
!===============================================================================
END MODULE shr_orb_mod