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

Last change on this file since 5296 was 5285, checked in by abarral, 3 days ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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