source: trunk/LMDZ.GENERIC/libf/phystd/stellarlong.F @ 537

Last change on this file since 537 was 253, checked in by emillour, 13 years ago

Generic GCM

  • Massive update to version 0.7

EM+RW

File size: 2.6 KB
RevLine 
[253]1      SUBROUTINE stellarlong(pday,pstellong)
2      IMPLICIT NONE
3
4c=======================================================================
5c
6c   Objet:
7c   ------
8c
9c      Calcul de la distance soleil-planete et de la declinaison
10c   en fonction du jour de l'annee.
11c
12c
13c   Methode:
14c   --------
15c
16c      Calcul complet de l'elipse
17c
18c   Interface:
19c   ----------
20c
21c      Uncommon comprenant les parametres orbitaux.
22c
23c   Arguments:
24c   ----------
25c
26c   Input:
27c   ------
28c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
29c
30c   Output:
31c   -------
32c   pdist_star     distance entre le soleil et la planete
33c                 ( en unite astronomique pour utiliser la constante
34c                  solaire terrestre 1370 Wm-2 )
35c   pdecli        declinaison ( en radians )
36c
37c=======================================================================
38c-----------------------------------------------------------------------
39c   Declarations:
40c   -------------
41
42#include "planete.h"
43#include "comcstfi.h"
44
45c arguments:
46c ----------
47
48      REAL pday,pdist_star,pdecli,pstellong
49      LOGICAL lwrite
50
51c Local:
52c ------
53
54      REAL zanom,xref,zx0,zdx,zteta,zz
55      INTEGER iter
56
57
58c-----------------------------------------------------------------------
59c calcul de l'angle polaire et de la distance au soleil :
60c -------------------------------------------------------
61
62c   Initialisation eventuelle:
63      if(.not.unitastr.gt.1.e-4) then
64         call iniorbit(249.22,206.66,669.,485.,25.2)
65      endif
66
67c  calcul de l'zanomalie moyenne
68
69      zz=(pday-peri_day)/year_day
70      pi=2.*asin(1.)
71      zanom=2.*pi*(zz-nint(zz))
72      xref=abs(zanom)
73
74c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
75c  methode de Newton
76
77      zx0=xref+e_elips*sin(xref)
78      DO 110 iter=1,10
79         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
80         if(abs(zdx).le.(1.e-7)) goto 120
81         zx0=zx0+zdx
82110   continue
83120   continue
84      zx0=zx0+zdx
85      if(zanom.lt.0.) zx0=-zx0
86
87c zteta est la longitude solaire
88
89      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
90
91      pstellong=zteta-timeperi
92
93      IF(pstellong.LT.0.) pstellong=pstellong+2.*pi
94      IF(pstellong.GT.2.*pi) pstellong=pstellong-2.*pi
95c-----------------------------------------------------------------------
96c   sorties eventuelles:
97c   ---------------------
98
99c     IF (lwrite) THEN
100c        PRINT*,'jour de l"annee   :',pday
101c        PRINT*,'distance au soleil (en unite astronomique) :',pdist_star
102c        PRINT*,'declinaison (en degres) :',pdecli*180./pi
103c     ENDIF
104
105      RETURN
106      END
Note: See TracBrowser for help on using the repository browser.