/[MITgcm]/MITgcm_contrib/darwin/pkg/radtrans/radtrans_sun2000.F
ViewVC logotype

Contents of /MITgcm_contrib/darwin/pkg/radtrans/radtrans_sun2000.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.1 - (show annotations) (download)
Wed Jun 9 15:59:38 2010 UTC (15 years, 1 month ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin_ckpt63j_20120217, ctrb_darwin_ckpt62w_20110426, ctrb_darwin_ckpt62k_20101005, ctrb_darwin_ckpt63m_20120506, ctrb_darwin_ckpt66l_20171025, ctrb_darwin_ckpt66a_20161020, ctrb_darwin_ckpt62l_20101006, ctrb_darwin_ckpt65n_20150729, ctrb_darwin_ckpt64v_20140411, ctrb_darwin_ckpt62q_20110124, ctrb_darwin_ckpt65b_20140812, ctrb_darwin_ckpt62m_20101024, ctrb_darwin_ckpt63s_20120908, ctrb_darwin_ckpt62z_20110622, ctrb_darwin_ckpt65y_20160801, ctrb_darwin_ckpt66o_20180209, ctrb_darwin_ckpt63q_20120731, ctrb_darwin_ckpt62s_20110217, ctrb_darwin_ckpt63k_20120317, ctrb_darwin_ckpt64y_20140622, ctrb_darwin_ckpt64b_20121224, ctrb_darwin_ckpt64g_20130503, ctrb_darwin_ckpt65c_20140830, ctrb_darwin_ckpt62o_20101125, ctrb_darwin_ckpt64r_20131210, ctrb_darwin_ckpt65x_20160612, ctrb_darwin_ckpt66j_20170815, ctrb_darwin_ckpt64j_20130704, ctrb_darwin_ckpt65u_20160315, ctrb_darwin_ckpt62i_20100729, ctrb_darwin_ckpt62t_20110305, ctrb_darwin_ckpt62h_20100707, ctrb_darwin_ckpt64t_20140202, ctrb_darwin_ckpt64_20121012, ctrb_darwin_ckpt64s_20140105, ctrb_darwin_ckpt66k_20171025, ctrb_darwin_ckpt63i_20120124, ctrb_darwin_ckpt64h_20130528, ctrb_darwin_ckpt64u_20140308, ctrb_darwin_ckpt64k_20130723, ctrb_darwin_ckpt63g_20111220, ctrb_darwin_ckpt63b_20110830, ctrb_darwin_ckpt65s_20160114, ctrb_darwin_ckpt64a_20121116, ctrb_darwin_ckpt62p_20101217, ctrb_darwin_ckpt63_20110728, ctrb_darwin_ckpt62n_20101109, ctrb_darwin_ckpt62r_20110124, ctrb_darwin_ckpt62v_20110412, ctrb_darwin_ckpt64i_20130622, ctrb_darwin_ckpt65h_20141217, ctrb_darwin_ckpt65a_20140728, ctrb_darwin_ckpt63o_20120629, ctrb_darwin_ckpt63a_20110804, ctrb_darwin_ckpt66i_20170718, ctrb_darwin_ckpt62u_20110405, ctrb_darwin_ckpt66n_20180118, ctrb_darwin_ckpt65l_20150504, ctrb_darwin_ckpt63c_20111011, ctrb_darwin_ckpt65o_20150914, ctrb_darwin_ckpt62v_20110406, ctrb_darwin_ckpt65f_20141014, ctrb_darwin_ckpt63p_20120707, ctrb_darwin_ckpt64q_20131118, ctrb_darwin_ckpt65i_20150123, ctrb_darwin_ckpt63d_20111107, ctrb_darwin_ckpt66b_20161219, ctrb_darwin_ckpt65t_20160221, ctrb_darwin_ckpt65k_20150402, ctrb_darwin_ckpt64p_20131024, ctrb_darwin_ckpt65e_20140929, ctrb_darwin_ckpt64c_20130120, ctrb_darwin_ckpt64q_20131024, ctrb_darwin_ckpt65m_20150615, ctrb_darwin_ckpt64m_20130820, ctrb_darwin_ckpt66f_20170407, ctrb_darwin_ckpt63r_20120817, ctrb_darwin_ckpt65w_20160512, ctrb_darwin_ckpt62x_20110513, ctrb_darwin_ckpt65d_20140915, ctrb_darwin_ckpt63f_20111201, ctrb_darwin_ckpt64o_20131024, ctrb_darwin_ckpt64f_20130405, ctrb_darwin_ckpt65r_20151221, ctrb_darwin_ckpt65q_20151118, ctrb_darwin_ckpt66g_20170424, ctrb_darwin_ckpt65j_20150225, ctrb_darwin_ckpt64w_20140502, ctrb_darwin_ckpt62j_20100901, ctrb_darwin_ckpt63e_20111107, ctrb_darwin_ckpt64d_20130219, ctrb_darwin_ckpt63l_20120405, ctrb_darwin_ckpt66c_20170121, ctrb_darwin_ckpt64x_20140524, ctrb_darwin_ckpt66h_20170602, ctrb_darwin_ckpt66d_20170214, ctrb_darwin_ckpt63h_20111230, ctrb_darwin_ckpt64z_20140711, ctrb_darwin_ckpt65v_20160409, ctrb_darwin_ckpt65p_20151023, ctrb_darwin_ckpt65g_20141120, ctrb_darwin_ckpt63n_20120604, ctrb_darwin_ckpt64e_20130305, ctrb_darwin_ckpt66m_20171213, ctrb_darwin_ckpt66e_20170314, ctrb_darwin_ckpt64l_20130806, ctrb_darwin_ckpt65z_20160929, ctrb_darwin_ckpt62y_20110526, ctrb_darwin_ckpt65_20140718, ctrb_darwin_ckpt64n_20130826, HEAD 
add radtrans package
1 C $Header$
2 C $Name$
3
4 #include "RADTRANS_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: RADTRANS_SUN2000
8
9 C !INTERFACE: ======================================================
10 subroutine radtrans_sun2000 (radeg, iyr, imon, iday, sec,
11 O sunvec, rs)
12
13 C !DESCRIPTION:
14 c This subroutine computes the Sun vector in geocentric inertial
15 c (equatorial) coodinates. It uses the model referenced in The
16 c Astronomical Almanac for 1984, Section S (Supplement) and documented
17 c in Exact closed-form geolocation algorithm for Earth survey
18 c sensors, by F.S. Patt and W.W. Gregg, Int. Journal of Remote
19 c Sensing, 1993. The accuracy of the Sun vector is approximately 0.1
20 c arcminute.
21 c
22 c Arguments:
23 c
24 c Name Type I/O Description
25 c --------------------------------------------------------
26 c IYR I*4 I Year, four digits (i.e, 1993)
27 c IDAY I*4 I Day of year (1-366)
28 c SEC R*8 I Seconds of day
29 c SUN(3) R*8 O Unit Sun vector in geocentric inertial
30 c coordinates of date
31 c RS R*8 O Magnitude of the Sun vector (AU)
32 c
33 c Subprograms referenced:
34 c
35 c JD Computes Julian day from calendar date
36 c EPHPARMS Computes mean solar longitude and anomaly and
37 c mean lunar lontitude and ascending node
38 c NUTATE Compute nutation corrections to lontitude and
39 c obliquity
40 c
41 c Coded by: Frederick S. Patt, GSC, November 2, 1992
42 c Modified to include Earth constants subroutine by W. Gregg,
43 c May 11, 1993.
44
45 C !USES: ===========================================================
46 IMPLICIT NONE
47 #include "RADTRANS_VARS.h"
48
49 C !INPUT PARAMETERS: ===============================================
50 INTEGER iyr, imon, iday
51 _RL radeg, sec
52 c INTEGER myThid
53
54 C !OUTPUT PARAMETERS: ==============================================
55 _RL sunvec(3), rs
56
57 C !FUNCTIONS: ======================================================
58 INTEGER radtrans_jd
59 EXTERNAL radtrans_jd
60
61 C !LOCAL VARIABLES: ================================================
62 INTEGER nt
63 _RL xk,rjd,t,xls,gs,xlm,omega,g2,g4,g5,dls,xlsg,xlsa
64 parameter (xk=0.0056932) !Constant of aberration
65 CEOP
66
67 c Compute floating point days since Jan 1.5, 2000
68 c Note that the Julian day starts at noon on the specified date
69 rjd = float(radtrans_jd(iyr,imon,iday))
70 t = rjd - 2451545.0D0 + (sec-43200.0D0)/86400.0D0
71
72 c Compute solar ephemeris parameters
73 call radtrans_ephparms (t, xls, gs, xlm, omega)
74
75 c Check if need to compute nutation corrections for this day
76 nt = int(t)
77 if (nt.ne.nutime) then
78 nutime = nt
79 call radtrans_nutate (radeg, t, xls, gs, xlm, omega, dpsi, eps)
80 end if
81
82 c Compute planet mean anomalies
83 c Venus Mean Anomaly
84 g2 = 50.40828D0 + 1.60213022D0*t
85 g2 = mod(g2,360.0)
86
87 c Mars Mean Anomaly
88 g4 = 19.38816D0 + 0.52402078D0*t
89 g4 = mod(g4,360.0)
90
91 c Jupiter Mean Anomaly
92 g5 = 20.35116D0 + 0.08309121D0*t
93 g5 = mod(g5,360.0)
94
95 c Compute solar distance (AU)
96 rs = 1.00014D0 - 0.01671D0*cos(gs/radeg)
97 * - 0.00014D0*cos(2.0D0*gs/radeg)
98
99 c Compute Geometric Solar Longitude
100 dls = (6893.0D0 - 4.6543463D-4*t)*sin(gs/radeg)
101 * + 72.0D0*sin(2.0D0*gs/radeg)
102 * - 7.0D0*cos((gs - g5)/radeg)
103 * + 6.0D0*sin((xlm - xls)/radeg)
104 * + 5.0D0*sin((4.0D0*gs - 8.0D0*g4 + 3.0D0*g5)/radeg)
105 * - 5.0D0*cos((2.0D0*gs - 2.0D0*g2)/radeg)
106 * - 4.0D0*sin((gs - g2)/radeg)
107 * + 4.0D0*cos((4.0D0*gs - 8.0D0*g4 + 3.0D0*g5)/radeg)
108 * + 3.0D0*sin((2.0D0*gs - 2.0D0*g2)/radeg)
109 * - 3.0D0*sin(g5/radeg)
110 * - 3.0D0*sin((2.0D0*gs - 2.0D0*g5)/radeg) !arcseconds
111
112 xlsg = xls + dls/3600.0D0
113
114 c Compute Apparent Solar Longitude; includes corrections for nutation
115 c in longitude and velocity aberration
116 xlsa = xlsg + dpsi - xk/rs
117
118 c Compute unit Sun vector
119 sunvec(1) = cos(xlsa/radeg)
120 sunvec(2) = sin(xlsa/radeg)*cos(eps/radeg)
121 sunvec(3) = sin(xlsa/radeg)*sin(eps/radeg)
122 c type *,' Sunlon = ',xlsg,xlsa,eps
123
124 return
125 end
126 c
  ViewVC Help
Powered by ViewVC 1.1.22