source: dynamico_lmdz/simple_physics/phyparam/param/solarlong.F @ 4182

Last change on this file since 4182 was 4176, checked in by dubos, 6 years ago

simple_physics : copy code from FH

File size: 2.4 KB
Line 
1      SUBROUTINE solarlong(pday,psollong)
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   lwrite        clef logique pour sorties de controle
30c
31c   Output:
32c   -------
33c   pdist_sol     distance entre le soleil et la planete
34c                 ( en unite astronomique pour utiliser la constante
35c                  solaire terrestre 1370 Wm-2 )
36c   pdecli        declinaison ( en radians )
37c
38c=======================================================================
39c-----------------------------------------------------------------------
40c   Declarations:
41c   -------------
42
43#include "comorbit.h"
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      zanom=2.*pi*(zz-nint(zz))
66      xref=abs(zanom)
67
68c  resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
69c  methode de Newton
70
71      zx0=xref+e_elips*sin(xref)
72      DO 110 iter=1,10
73         zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
74         if(abs(zdx).le.(1.e-7)) goto 120
75         zx0=zx0+zdx
76110   continue
77120   continue
78      zx0=zx0+zdx
79      if(zanom.lt.0.) zx0=-zx0
80
81c zteta est la longitude solaire
82
83      zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
84
85      psollong=zteta-timeperi
86
87      IF(psollong.LT.0.) psollong=psollong+2.*pi
88      IF(psollong.GT.2.*pi) psollong=psollong-2.*pi
89c-----------------------------------------------------------------------
90c   sorties eventuelles:
91c   ---------------------
92
93      IF (lwrite) THEN
94         PRINT*,'jour de l"annee   :',pday
95         PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
96         PRINT*,'declinaison (en degres) :',pdecli*180./pi
97      ENDIF
98
99      RETURN
100      END
Note: See TracBrowser for help on using the repository browser.