source: LMDZ6/branches/Amaury_dev/libf/phylmd/solarlong.F90 @ 5142

Last change on this file since 5142 was 5142, checked in by abarral, 3 months ago

Put cvparam.h, fcg_gcssold.h, planete.h, tsoilnudge.h, YOECUMF.h into 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.5 KB
Line 
1SUBROUTINE solarlong(pday, psollong, pdist_sol)
2
3  USE ioipsl
4  USE lmdz_print_control, ONLY: lunout
5  USE lmdz_planete, ONLY: aphelie, periheli, year_day, peri_day, obliquit, timeperi, e_elips, p_elips, unitastr
6
7  IMPLICIT 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  include "YOMCST.h"
48
49  ! arguments:
50  ! ----------
51
52  REAL pday, pdist_sol, pdecli, psollong
53  LOGICAL lwrite
54
55  ! Local:
56  ! ------
57
58  REAL zanom, xref, zx0, zdx, zteta, zz, pi
59  INTEGER iter
60  REAL :: pyear_day, pperi_day
61  REAL :: jd_eq, jd_peri
62  LOGICAL, SAVE :: first = .TRUE.
63  !$OMP THREADPRIVATE(first)
64
65  ! -----------------------------------------------------------------------
66  ! calcul de l'angle polaire et de la distance au soleil :
67  ! -------------------------------------------------------
68
69  ! Initialisation eventuelle:
70  IF (first) THEN
71    CALL ioget_calendar(pyear_day)
72    CALL ymds2ju(2000, 3, 21, 0., jd_eq)
73    CALL ymds2ju(2001, 1, 4, 0., jd_peri)
74    pperi_day = jd_peri - jd_eq
75    pperi_day = r_peri + 180.
76    WRITE (lunout, *) ' Number of days in a year = ', pyear_day
77    ! CALL iniorbit(249.22,206.66,669.,485.,25.2)
78    CALL iniorbit(152.59, 146.61, pyear_day, pperi_day, r_incl)
79    first = .FALSE.
80  END IF
81
82  ! calcul de l'zanomalie moyenne
83
84  zz = (pday - peri_day) / year_day
85  pi = 2. * asin(1.)
86  zanom = 2. * pi * (zz - nint(zz))
87  xref = abs(zanom)
88
89  ! resolution de l'equation horaire  zx0 - e * sin (zx0) = xref
90  ! methode de Newton
91
92  ! zx0=xref+e_elips*sin(xref)
93  zx0 = xref + r_ecc * sin(xref)
94  DO iter = 1, 10
95    ! zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0))
96    zdx = -(zx0 - r_ecc * sin(zx0) - xref) / (1. - r_ecc * cos(zx0))
97    IF (abs(zdx)<=(1.E-7)) GO TO 120
98    zx0 = zx0 + zdx
99  END DO
100  120 CONTINUE
101  zx0 = zx0 + zdx
102  IF (zanom<0.) zx0 = -zx0
103
104  ! zteta est la longitude solaire
105
106  ! zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
107  zteta = 2. * atan(sqrt((1. + r_ecc) / (1. - r_ecc)) * tan(zx0 / 2.))
108
109  psollong = zteta - timeperi
110
111  IF (psollong<0.) psollong = psollong + 2. * pi
112  IF (psollong>2. * pi) psollong = psollong - 2. * pi
113
114  psollong = psollong * 180. / pi
115
116  ! distance soleil
117
118  pdist_sol = (1 - r_ecc * r_ecc) / (1 + r_ecc * cos(pi / 180. * (psollong - &
119          (r_peri + 180.0))))
120  ! pdist_sol = (1-e_elips*e_elips)
121  ! &      /(1+e_elips*COS(pi/180.*(psollong-(R_peri+180.0))))
122  ! -----------------------------------------------------------------------
123  ! sorties eventuelles:
124  ! ---------------------
125
126  ! IF (lwrite) THEN
127  ! PRINT*,'jour de l"annee   :',pday
128  ! PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol
129  ! PRINT*,'declinaison (en degres) :',pdecli*180./pi
130  ! ENDIF
131
132END SUBROUTINE solarlong
Note: See TracBrowser for help on using the repository browser.