source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/iniorbit.F @ 4999

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

Modifications nécessaires a l'inclusion d'un calendrier réaliste.
La date courante est calculée dans leapfrog.F et exprimée en Jour Julien
(modifié). On en a profité pour faire un peu de ménage dans la gestion des dates
du modèle.
Dans la physique, on utilise les routines de passages entre calendrier Julien et
Gregorien incluses dans IOIPSL pour calculer le nombre de jours écoulés depuis le
1er janvier (pour les conditions aux limites) ou l'equinoxe (pour le calcul de
la longitude solaire). Le calcul de l'orbite reprend celui du gcm planétaire
(codé par FH)
On décide du calendrier à utiliser à l'aide du paramètre calend du run.def. Par
défaut celui-ci est à earth_360d
LF

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.