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

Last change on this file since 2537 was 1384, checked in by emillour, 10 years ago

Generic GCM:

  • Some code cleanup: turning comcstfi.h into module comcstfi_mod.F90

EM

File size: 2.5 KB
Line 
1      SUBROUTINE stellarlong(pday,pstellong)
2     
3      USE planete_mod, ONLY: year_day, peri_day, e_elips, timeperi
4      use comcstfi_mod, only: pi
5      IMPLICIT NONE
6
7c=======================================================================
8c
9c   Objet:
10c   ------
11c
12c      Calcul de la distance soleil-planete et de la declinaison
13c   en fonction du jour de l'annee.
14c
15c
16c   Methode:
17c   --------
18c
19c      Calcul complet de l'elipse
20c
21c   Interface:
22c   ----------
23c
24c      Uncommon comprenant les parametres orbitaux.
25c
26c   Arguments:
27c   ----------
28c
29c   Input:
30c   ------
31c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
32c
33c   Output:
34c   -------
35c   pdist_star     distance entre le soleil et la planete
36c                 ( en unite astronomique pour utiliser la constante
37c                  solaire terrestre 1370 Wm-2 )
38c   pdecli        declinaison ( en radians )
39c
40c=======================================================================
41c-----------------------------------------------------------------------
42c   Declarations:
43c   -------------
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  calcul de l'zanomalie moyenne
63
64      zz=(pday-peri_day)/year_day
65      pi=2.*asin(1.)
66      zanom=2.*pi*(zz-nint(zz))
67      xref=abs(zanom)
68
69c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
70c  methode de Newton
71
72      zx0=xref+e_elips*sin(xref)
73      DO 110 iter=1,10
74         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
75         if(abs(zdx).le.(1.e-7)) goto 120
76         zx0=zx0+zdx
77110   continue
78120   continue
79      zx0=zx0+zdx
80      if(zanom.lt.0.) zx0=-zx0
81
82c zteta est la longitude solaire
83
84      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
85
86      pstellong=zteta-timeperi
87
88      IF(pstellong.LT.0.) pstellong=pstellong+2.*pi
89      IF(pstellong.GT.2.*pi) pstellong=pstellong-2.*pi
90c-----------------------------------------------------------------------
91c   sorties eventuelles:
92c   ---------------------
93
94c     IF (lwrite) THEN
95c        PRINT*,'jour de l"annee   :',pday
96c        PRINT*,'distance au soleil (en unite astronomique) :',pdist_star
97c        PRINT*,'declinaison (en degres) :',pdecli*180./pi
98c     ENDIF
99
100      RETURN
101      END
Note: See TracBrowser for help on using the repository browser.