source: LMDZ6/trunk/libf/phylmd/solarlong.f90 @ 5296

Last change on this file since 5296 was 5285, checked in by abarral, 3 days ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules

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