source: LMDZ5/branches/testing/libf/phylmd/iniorbit.F90 @ 3525

Last change on this file since 3525 was 1999, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1920:1997 into testing branch

  • 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  IMPLICIT NONE
3
4  ! =======================================================================
5
6  ! Auteur:
7  ! -------
8  ! Frederic Hourdin      22 Fevrier 1991
9
10  ! Objet:
11  ! ------
12  ! Initialisation du sous programme orbite qui calcule
13  ! a une date donnee de l'annee de duree year_day commencant
14  ! a l'equinoxe de printemps et dont le perihelie se situe
15  ! a la date peri_day, la distance au soleil et la declinaison.
16
17  ! Interface:
18  ! ----------
19  ! - Doit etre appele avant d'utiliser orbite.
20  ! - initialise une partie du common planete.h
21
22  ! Arguments:
23  ! ----------
24
25  ! Input:
26  ! ------
27  ! aphelie       \   aphelie et perihelie de l'orbite
28  ! periheli      /   en millions de kilometres.
29
30  ! =======================================================================
31
32  ! -----------------------------------------------------------------------
33  ! Declarations:
34  ! -------------
35
36  include "planete.h"
37  include "YOMCST.h"
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.