source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/iniorbit.F @ 1373

Last change on this file since 1373 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

File size: 2.7 KB
Line 
1      SUBROUTINE iniorbit
2     $     (paphelie,pperiheli,pyear_day,pperi_day,pobliq)
3      IMPLICIT NONE
4
5c=======================================================================
6c
7c   Auteur:
8c   -------
9c     Frederic Hourdin      22 Fevrier 1991
10c
11c   Objet:
12c   ------
13c    Initialisation du sous programme orbite qui calcule
14c    a une date donnee de l'annee de duree year_day commencant
15c    a l'equinoxe de printemps et dont le perihelie se situe
16c    a la date peri_day, la distance au soleil et la declinaison.
17c
18c   Interface:
19c   ----------
20c   - Doit etre appele avant d'utiliser orbite.
21c   - initialise une partie du common planete.h
22c
23c   Arguments:
24c   ----------
25c
26c   Input:
27c   ------
28c   aphelie       \   aphelie et perihelie de l'orbite
29c   periheli      /   en millions de kilometres.
30c
31c=======================================================================
32
33c-----------------------------------------------------------------------
34c   Declarations:
35c   -------------
36
37#include "planete.h"
38#include "YOMCST.h"
39
40c   Arguments:
41c   ----------
42
43      REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq
44
45c   Local:
46c   ------
47
48      REAL zxref,zanom,zz,zx0,zdx, pi
49      INTEGER iter
50
51c-----------------------------------------------------------------------
52
53      pi=2.*asin(1.)
54
55      aphelie =paphelie
56      periheli=pperiheli
57      year_day=pyear_day
58      obliquit=pobliq
59      peri_day=pperi_day
60
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
69
70      print*,'e_elips',e_elips
71      print*,'p_elips',p_elips
72
73c-----------------------------------------------------------------------
74c calcul de l'angle polaire et de la distance au soleil :
75c -------------------------------------------------------
76
77c  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
84c  resolution de l'equation horaire  zx0 - e * sin (zx0) = zxref
85c  methode de Newton
86
87      zx0=zxref+R_ecc*sin(zxref)
88      DO 110 iter=1,100
89         zdx=-(zx0-R_ecc*sin(zx0)-zxref)/(1.-R_ecc*cos(zx0))
90         if(abs(zdx).le.(1.e-12)) goto 120
91         zx0=zx0+zdx
92110   continue
93120   continue
94      zx0=zx0+zdx
95      if(zanom.lt.0.) zx0=-zx0
96      PRINT*,'zx0   ',zx0
97
98c 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
104      END
Note: See TracBrowser for help on using the repository browser.