source: LMDZ6/trunk/libf/phylmd/iniorbit.f90 @ 5407

Last change on this file since 5407 was 5300, checked in by abarral, 6 weeks ago

Fix remaining planete.h

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 2.6 KB
Line 
1SUBROUTINE iniorbit(paphelie, pperiheli, pyear_day, pperi_day, pobliq)
2  USE yomcst_mod_h
3  USE planete_mod_h
4IMPLICIT NONE
5
6  ! =======================================================================
7
8  ! Auteur:
9  ! -------
10  ! Frederic Hourdin      22 Fevrier 1991
11
12  ! Objet:
13  ! ------
14  ! Initialisation du sous programme orbite qui calcule
15  ! a une date donnee de l'annee de duree year_day commencant
16  ! a l'equinoxe de printemps et dont le perihelie se situe
17  ! a la date peri_day, la distance au soleil et la declinaison.
18
19  ! Interface:
20  ! ----------
21  ! - Doit etre appele avant d'utiliser orbite.
22  ! - initialise une partie du common planete.h
23
24  ! Arguments:
25  ! ----------
26
27  ! Input:
28  ! ------
29  ! aphelie       \   aphelie et perihelie de l'orbite
30  ! periheli      /   en millions de kilometres.
31
32  ! =======================================================================
33
34  ! -----------------------------------------------------------------------
35  ! Declarations:
36  ! -------------
37
38
39  ! Arguments:
40  ! ----------
41
42  REAL paphelie, pperiheli, pyear_day, pperi_day, pobliq
43
44  ! Local:
45  ! ------
46
47  REAL zxref, zanom, zz, zx0, zdx, pi
48  INTEGER iter
49
50  ! -----------------------------------------------------------------------
51
52  pi = 2.*asin(1.)
53
54  aphelie = paphelie
55  periheli = pperiheli
56  year_day = pyear_day
57  obliquit = pobliq
58  peri_day = pperi_day
59
60  PRINT *, 'Perihelie en Mkm  ', periheli
61  PRINT *, 'Aphelie  en Mkm   ', aphelie
62  PRINT *, 'obliquite en degres  :', obliquit
63  PRINT *, 'Jours dans l annee : ', year_day
64  PRINT *, 'Date perihelie : ', peri_day
65  unitastr = 149.597870
66  e_elips = (aphelie-periheli)/(periheli+aphelie)
67  p_elips = 0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr
68
69  PRINT *, 'e_elips', e_elips
70  PRINT *, 'p_elips', p_elips
71
72  ! -----------------------------------------------------------------------
73  ! calcul de l'angle polaire et de la distance au soleil :
74  ! -------------------------------------------------------
75
76  ! calcul de l'zanomalie moyenne
77
78  zz = (year_day-pperi_day)/year_day
79  zanom = 2.*pi*(zz-nint(zz))
80  zxref = abs(zanom)
81  PRINT *, 'zanom  ', zanom
82
83  ! resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
84  ! methode de Newton
85
86  zx0 = zxref + r_ecc*sin(zxref)
87  DO iter = 1, 100
88    zdx = -(zx0-r_ecc*sin(zx0)-zxref)/(1.-r_ecc*cos(zx0))
89    IF (abs(zdx)<=(1.E-12)) GO TO 120
90    zx0 = zx0 + zdx
91  END DO
92120 CONTINUE
93  zx0 = zx0 + zdx
94  IF (zanom<0.) zx0 = -zx0
95  PRINT *, 'zx0   ', zx0
96
97  ! zteta est la longitude solaire
98
99  timeperi = 2.*atan(sqrt((1.+r_ecc)/(1.-r_ecc))*tan(zx0/2.))
100  PRINT *, 'longitude solaire du perihelie timeperi = ', timeperi
101
102  RETURN
103END SUBROUTINE iniorbit
Note: See TracBrowser for help on using the repository browser.