source: trunk/LMDZ.MARS/libf/phymars/solarlong.F @ 2616

Last change on this file since 2616 was 1226, checked in by aslmd, 11 years ago

LMDZ.MARS : Replaced comcstfi and planete includes by modules.

File size: 2.5 KB
Line 
1      SUBROUTINE solarlong(pday,psollong)
2      use planete_h
3      USE comcstfi_h
4      IMPLICIT NONE
5
6c=======================================================================
7c
8c   Objet:
9c   ------
10c
11c      Calcul de la distance soleil-planete et de la declinaison
12c   en fonction du jour de l'annee.
13c
14c
15c   Methode:
16c   --------
17c
18c      Calcul complet de l'elipse
19c
20c   Interface:
21c   ----------
22c
23c      Uncommon comprenant les parametres orbitaux.
24c
25c   Arguments:
26c   ----------
27c
28c   Input:
29c   ------
30c   pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
31c   lwrite        clef logique pour sorties de controle
32c
33c   Output:
34c   -------
35c   pdist_sol     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_sol,pdecli,psollong
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      psollong=zteta-timeperi
87
88      IF(psollong.LT.0.) psollong=psollong+2.*pi
89      IF(psollong.GT.2.*pi) psollong=psollong-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_sol
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.