source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/solarlong.F90 @ 3069

Last change on this file since 3069 was 1999, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1920:1997 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 3.3 KB
Line 
1SUBROUTINE solarlong(pday, psollong, pdist_sol)
2
3  USE ioipsl
4
5  IMPLICIT NONE
6
7  ! =======================================================================
8
9  ! Objet:
10  ! ------
11
12  ! Calcul de la distance soleil-planete et de la declinaison
13  ! en fonction du jour de l'annee.
14
15
16  ! Methode:
17  ! --------
18
19  ! Calcul complet de l'elipse
20
21  ! Interface:
22  ! ----------
23
24  ! Uncommon comprenant les parametres orbitaux.
25
26  ! Arguments:
27  ! ----------
28
29  ! Input:
30  ! ------
31  ! pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
32  ! lwrite        clef logique pour sorties de controle
33
34  ! Output:
35  ! -------
36  ! pdist_sol     distance entre le soleil et la planete
37  ! ( en unite astronomique pour utiliser la constante
38  ! solaire terrestre 1370 Wm-2 )
39  ! pdecli        declinaison ( en radians )
40
41  ! =======================================================================
42  ! -----------------------------------------------------------------------
43  ! Declarations:
44  ! -------------
45
46  include "planete.h"
47  include "YOMCST.h"
48  include 'iniprint.h'
49
50  ! arguments:
51  ! ----------
52
53  REAL pday, pdist_sol, pdecli, psollong
54  LOGICAL lwrite
55
56  ! Local:
57  ! ------
58
59  REAL zanom, xref, zx0, zdx, zteta, zz, pi
60  INTEGER iter
61  REAL :: pyear_day, pperi_day
62  REAL :: jd_eq, jd_peri
63  LOGICAL, SAVE :: first = .TRUE.
64  !$OMP THREADPRIVATE(first)
65
66  ! -----------------------------------------------------------------------
67  ! calcul de l'angle polaire et de la distance au soleil :
68  ! -------------------------------------------------------
69
70  ! Initialisation eventuelle:
71  IF (first) THEN
72    CALL ioget_calendar(pyear_day)
73    CALL ymds2ju(2000, 3, 21, 0., jd_eq)
74    CALL ymds2ju(2001, 1, 4, 0., jd_peri)
75    pperi_day = jd_peri - jd_eq
76    pperi_day = r_peri + 180.
77    WRITE (lunout, *) ' Number of days in a year = ', pyear_day
78    ! call iniorbit(249.22,206.66,669.,485.,25.2)
79    CALL iniorbit(152.59, 146.61, pyear_day, pperi_day, r_incl)
80    first = .FALSE.
81  END IF
82
83  ! calcul de l'zanomalie moyenne
84
85  zz = (pday-peri_day)/year_day
86  pi = 2.*asin(1.)
87  zanom = 2.*pi*(zz-nint(zz))
88  xref = abs(zanom)
89
90  ! resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
91  ! methode de Newton
92
93  ! zx0=xref+e_elips*sin(xref)
94  zx0 = xref + r_ecc*sin(xref)
95  DO iter = 1, 10
96    ! zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
97    zdx = -(zx0-r_ecc*sin(zx0)-xref)/(1.-r_ecc*cos(zx0))
98    IF (abs(zdx)<=(1.E-7)) GO TO 120
99    zx0 = zx0 + zdx
100  END DO
101120 CONTINUE
102  zx0 = zx0 + zdx
103  IF (zanom<0.) zx0 = -zx0
104
105  ! zteta est la longitude solaire
106
107  ! zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
108  zteta = 2.*atan(sqrt((1.+r_ecc)/(1.-r_ecc))*tan(zx0/2.))
109
110  psollong = zteta - timeperi
111
112  IF (psollong<0.) psollong = psollong + 2.*pi
113  IF (psollong>2.*pi) psollong = psollong - 2.*pi
114
115  psollong = psollong*180./pi
116
117  ! distance soleil
118
119  pdist_sol = (1-r_ecc*r_ecc)/(1+r_ecc*cos(pi/180.*(psollong- &
120    (r_peri+180.0))))
121  ! pdist_sol = (1-e_elips*e_elips)
122  ! &      /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0))))
123  ! -----------------------------------------------------------------------
124  ! sorties eventuelles:
125  ! ---------------------
126
127  ! IF (lwrite) THEN
128  ! PRINT*,'jour de l"annee   :',pday
129  ! PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
130  ! PRINT*,'declinaison (en degres) :',pdecli*180./pi
131  ! ENDIF
132
133  RETURN
134END SUBROUTINE solarlong
Note: See TracBrowser for help on using the repository browser.