source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/solarlong.F90 @ 5376

Last change on this file since 5376 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM

  • 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  IMPLICIT NONE
7
8  ! =======================================================================
9
10  ! Objet:
11  ! ------
12
13  ! Calcul de la distance soleil-planete et de la declinaison
14  ! en fonction du jour de l'annee.
15
16
17  ! Methode:
18  ! --------
19
20  ! Calcul complet de l'elipse
21
22  ! Interface:
23  ! ----------
24
25  ! Uncommon comprenant les parametres orbitaux.
26
27  ! Arguments:
28  ! ----------
29
30  ! Input:
31  ! ------
32  ! pday          jour de l'annee (le jour 0 correspondant a l'equinoxe)
33  ! lwrite        clef logique pour sorties de controle
34
35  ! Output:
36  ! -------
37  ! pdist_sol     distance entre le soleil et la planete
38  ! ( en unite astronomique pour utiliser la constante
39  ! solaire terrestre 1370 Wm-2 )
40  ! pdecli        declinaison ( en radians )
41
42  ! =======================================================================
43  ! -----------------------------------------------------------------------
44  ! Declarations:
45  ! -------------
46
47  include "planete.h"
48  include "YOMCST.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.