source: LMDZ4/branches/LMDZ4V5.0-LF/libf/phylmd/solarlong.F @ 5394

Last change on this file since 5394 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

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