1 |
|
1153 |
SUBROUTINE solarlong(pday, psollong, pdist_sol) |
2 |
|
|
|
3 |
|
|
USE ioipsl |
4 |
|
|
USE print_control_mod, ONLY: lunout |
5 |
|
|
|
6 |
|
|
IMPLICIT NONE |
7 |
|
|
|
8 |
|
|
! ======================================================================= |
9 |
|
|
|
10 |
|
|
! Objet: |
11 |
|
|
! ------ |
12 |
|
|
|
13 |
|
|
! Calcul de la distance soleil-planete et de la declinaison |
14 |
|
|
! en fonction du jour de l'annee. |
15 |
|
|
|
16 |
|
|
|
17 |
|
|
! Methode: |
18 |
|
|
! -------- |
19 |
|
|
|
20 |
|
|
! Calcul complet de l'elipse |
21 |
|
|
|
22 |
|
|
! Interface: |
23 |
|
|
! ---------- |
24 |
|
|
|
25 |
|
|
! Uncommon comprenant les parametres orbitaux. |
26 |
|
|
|
27 |
|
|
! Arguments: |
28 |
|
|
! ---------- |
29 |
|
|
|
30 |
|
|
! Input: |
31 |
|
|
! ------ |
32 |
|
|
! pday jour de l'annee (le jour 0 correspondant a l'equinoxe) |
33 |
|
|
! lwrite clef logique pour sorties de controle |
34 |
|
|
|
35 |
|
|
! Output: |
36 |
|
|
! ------- |
37 |
|
|
! pdist_sol distance entre le soleil et la planete |
38 |
|
|
! ( en unite astronomique pour utiliser la constante |
39 |
|
|
! solaire terrestre 1370 Wm-2 ) |
40 |
|
|
! pdecli declinaison ( en radians ) |
41 |
|
|
|
42 |
|
|
! ======================================================================= |
43 |
|
|
! ----------------------------------------------------------------------- |
44 |
|
|
! Declarations: |
45 |
|
|
! ------------- |
46 |
|
|
|
47 |
|
|
include "planete.h" |
48 |
|
|
include "YOMCST.h" |
49 |
|
|
|
50 |
|
|
! arguments: |
51 |
|
|
! ---------- |
52 |
|
|
|
53 |
|
|
REAL pday, pdist_sol, pdecli, psollong |
54 |
|
|
LOGICAL lwrite |
55 |
|
|
|
56 |
|
|
! Local: |
57 |
|
|
! ------ |
58 |
|
|
|
59 |
|
|
REAL zanom, xref, zx0, zdx, zteta, zz, pi |
60 |
|
|
INTEGER iter |
61 |
|
|
REAL :: pyear_day, pperi_day |
62 |
|
|
REAL :: jd_eq, jd_peri |
63 |
|
|
LOGICAL, SAVE :: first = .TRUE. |
64 |
|
|
!$OMP THREADPRIVATE(first) |
65 |
|
|
|
66 |
|
|
! ----------------------------------------------------------------------- |
67 |
|
|
! calcul de l'angle polaire et de la distance au soleil : |
68 |
|
|
! ------------------------------------------------------- |
69 |
|
|
|
70 |
|
|
! Initialisation eventuelle: |
71 |
✓✓ |
288 |
IF (first) THEN |
72 |
|
1 |
CALL ioget_calendar(pyear_day) |
73 |
|
1 |
CALL ymds2ju(2000, 3, 21, 0., jd_eq) |
74 |
|
1 |
CALL ymds2ju(2001, 1, 4, 0., jd_peri) |
75 |
|
|
pperi_day = jd_peri - jd_eq |
76 |
|
1 |
pperi_day = r_peri + 180. |
77 |
|
1 |
WRITE (lunout, *) ' Number of days in a year = ', pyear_day |
78 |
|
|
! call iniorbit(249.22,206.66,669.,485.,25.2) |
79 |
|
1 |
CALL iniorbit(152.59, 146.61, pyear_day, pperi_day, r_incl) |
80 |
|
1 |
first = .FALSE. |
81 |
|
|
END IF |
82 |
|
|
|
83 |
|
|
! calcul de l'zanomalie moyenne |
84 |
|
|
|
85 |
|
288 |
zz = (pday-peri_day)/year_day |
86 |
|
|
pi = 2.*asin(1.) |
87 |
|
288 |
zanom = 2.*pi*(zz-nint(zz)) |
88 |
|
288 |
xref = abs(zanom) |
89 |
|
|
|
90 |
|
|
! resolution de l'equation horaire zx0 - e * sin (zx0) = xref |
91 |
|
|
! methode de Newton |
92 |
|
|
|
93 |
|
|
! zx0=xref+e_elips*sin(xref) |
94 |
|
288 |
zx0 = xref + r_ecc*sin(xref) |
95 |
✓✗ |
572 |
DO iter = 1, 10 |
96 |
|
|
! zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) |
97 |
|
572 |
zdx = -(zx0-r_ecc*sin(zx0)-xref)/(1.-r_ecc*cos(zx0)) |
98 |
✓✓ |
572 |
IF (abs(zdx)<=(1.E-7)) GO TO 120 |
99 |
|
288 |
zx0 = zx0 + zdx |
100 |
|
|
END DO |
101 |
|
|
120 CONTINUE |
102 |
|
288 |
zx0 = zx0 + zdx |
103 |
✓✓ |
288 |
IF (zanom<0.) zx0 = -zx0 |
104 |
|
|
|
105 |
|
|
! zteta est la longitude solaire |
106 |
|
|
|
107 |
|
|
! zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) |
108 |
|
288 |
zteta = 2.*atan(sqrt((1.+r_ecc)/(1.-r_ecc))*tan(zx0/2.)) |
109 |
|
|
|
110 |
|
288 |
psollong = zteta - timeperi |
111 |
|
|
|
112 |
✓✗ |
288 |
IF (psollong<0.) psollong = psollong + 2.*pi |
113 |
✗✓ |
288 |
IF (psollong>2.*pi) psollong = psollong - 2.*pi |
114 |
|
|
|
115 |
|
288 |
psollong = psollong*180./pi |
116 |
|
|
|
117 |
|
|
! distance soleil |
118 |
|
|
|
119 |
|
|
pdist_sol = (1-r_ecc*r_ecc)/(1+r_ecc*cos(pi/180.*(psollong- & |
120 |
|
288 |
(r_peri+180.0)))) |
121 |
|
|
! pdist_sol = (1-e_elips*e_elips) |
122 |
|
|
! & /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0)))) |
123 |
|
|
! ----------------------------------------------------------------------- |
124 |
|
|
! sorties eventuelles: |
125 |
|
|
! --------------------- |
126 |
|
|
|
127 |
|
|
! IF (lwrite) THEN |
128 |
|
|
! PRINT*,'jour de l"annee :',pday |
129 |
|
|
! PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol |
130 |
|
|
! PRINT*,'declinaison (en degres) :',pdecli*180./pi |
131 |
|
|
! ENDIF |
132 |
|
|
|
133 |
|
288 |
RETURN |
134 |
|
|
END SUBROUTINE solarlong |