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