My Project
 All Classes Files Functions Variables Macros
orbite.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4 c======================================================================
5  SUBROUTINE orbite(xjour,longi,dist)
6  IMPLICIT none
7 c======================================================================
8 c Auteur(s): Z.X. Li (LMD/CNRS) (adapte du GCM du LMD) date: 19930818
9 c Objet: pour un jour donne, calculer la longitude vraie de la terre
10 c (par rapport au point vernal-21 mars) dans son orbite solaire
11 c calculer aussi la distance terre-soleil (unite astronomique)
12 c======================================================================
13 c Arguments:
14 c xjour--INPUT--R- jour de l'annee a compter du 1er janvier
15 c longi--OUTPUT-R- longitude vraie en degres par rapport au point
16 c vernal (21 mars) en degres
17 c dist---OUTPUT-R- distance terre-soleil (par rapport a la moyenne)
18  REAL xjour, longi, dist
19 c======================================================================
20 #include "YOMCST.h"
21 C
22 C -- Variables dynamiques locales
23  REAL pir,xl,xllp,xee,xse,xlam,dlamm,anm,ranm,anv,ranv
24 C
25  pir = 4.0*atan(1.0) / 180.0
26  xl=r_peri+180.0
27  xllp=xl*pir
28  xee=r_ecc*r_ecc
29  xse=sqrt(1.0-xee)
30  xlam = (r_ecc/2.0+r_ecc*xee/8.0)*(1.0+xse)*sin(xllp)
31  . - xee/4.0*(0.5+xse)*sin(2.0*xllp)
32  . + r_ecc*xee/8.0*(1.0/3.0+xse)*sin(3.0*xllp)
33  xlam=2.0*xlam/pir
34  dlamm=xlam+(xjour-81.0)
35  anm=dlamm-xl
36  ranm=anm*pir
37  xee=xee*r_ecc
38  ranv=ranm+(2.0*r_ecc-xee/4.0)*sin(ranm)
39  . +5.0/4.0*r_ecc*r_ecc*sin(2.0*ranm)
40  . +13.0/12.0*xee*sin(3.0*ranm)
41 c
42  anv=ranv/pir
43  longi=anv+xl
44 C
45  dist = (1-r_ecc*r_ecc)
46  . /(1+r_ecc*cos(pir*(longi-(r_peri+180.0))))
47  RETURN
48  END
49 c======================================================================
50  SUBROUTINE angle(longi, lati, frac, muzero)
51  USE dimphy
52  IMPLICIT none
53 c======================================================================
54 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
55 c Objet: Calculer la duree d'ensoleillement pour un jour et la hauteur
56 c du soleil (cosinus de l'angle zinithal) moyenne sur la journee
57 c======================================================================
58 c Arguments:
59 c longi----INPUT-R- la longitude vraie de la terre dans son plan
60 c solaire a partir de l'equinoxe de printemps (degre)
61 c lati-----INPUT-R- la latitude d'un point sur la terre (degre)
62 c frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee
63 c par 24 heures (unite en fraction de 0 a 1)
64 c muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur
65 c la journee (0 a 1)
66 c======================================================================
67 cym#include "dimensions.h"
68 cym#include "dimphy.h"
69  REAL longi
70  REAL lati(klon), frac(klon), muzero(klon)
71 #include "YOMCST.h"
72  REAL lat, omega, lon_sun, lat_sun
73  REAL pi_local, incl
74  INTEGER i
75 c
76  pi_local = 4.0 * atan(1.0)
77  incl=r_incl * pi_local / 180.
78 c
79  lon_sun = longi * pi_local / 180.0
80  lat_sun = asin(sin(lon_sun)*sin(incl) )
81 c
82  DO i = 1, klon
83  lat = lati(i) * pi_local / 180.0
84 c
85  IF ( lat .GE. (pi_local/2.+lat_sun)
86  . .OR. lat.LE.(-pi_local/2.+lat_sun)) THEN
87  omega = 0.0 ! nuit polaire
88  ELSE IF ( lat.GE.(pi_local/2.-lat_sun)
89  . .OR. lat.LE.(-pi_local/2.-lat_sun)) THEN
90  omega = pi_local ! journee polaire
91  ELSE
92  omega = -tan(lat)*tan(lat_sun)
93  omega = acos(omega)
94  ENDIF
95 c
96  frac(i) = omega / pi_local
97 c
98  IF (omega .GT. 0.0) THEN
99  muzero(i) = sin(lat)*sin(lat_sun)
100  . + cos(lat)*cos(lat_sun)*sin(omega) / omega
101  ELSE
102  muzero(i) = 0.0
103  ENDIF
104  ENDDO
105 c
106  RETURN
107  END
108 c====================================================================
109  SUBROUTINE zenang(longi,gmtime,pdtrad,lat,long,
110  s pmu0,frac)
111  USE dimphy
112  IMPLICIT none
113 c=============================================================
114 c Auteur : O. Boucher (LMD/CNRS)
115 c d'apres les routines zenith et angle de Z.X. Li
116 c Objet : calculer les valeurs moyennes du cos de l'angle zenithal
117 c et l'ensoleillement moyen entre gmtime1 et gmtime2
118 c connaissant la declinaison, la latitude et la longitude.
119 c Rque : Different de la routine angle en ce sens que zenang
120 c fournit des moyennes de pmu0 et non des valeurs
121 c instantanees, du coup frac prend toutes les valeurs
122 c entre 0 et 1.
123 c Date : premiere version le 13 decembre 1994
124 c revu pour GCM le 30 septembre 1996
125 c===============================================================
126 c longi : la longitude vraie de la terre dans son plan
127 c solaire a partir de l'equinoxe de printemps (degre)
128 c gmtime : temps universel en fraction de jour
129 c pdtrad : pas de temps du rayonnement (secondes)
130 c lat------INPUT : latitude en degres
131 c long-----INPUT : longitude en degres
132 c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
133 c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
134 c================================================================
135 cym#include "dimensions.h"
136 cym#include "dimphy.h"
137 #include "YOMCST.h"
138 c================================================================
139  real, intent(in):: longi, gmtime, pdtrad
140  real lat(klon), long(klon), pmu0(klon), frac(klon)
141 c================================================================
142  integer i
143  real gmtime1, gmtime2
144  real pi_local, deux_pi_local, incl
145  real omega1, omega2, omega
146 c omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.
147 c omega : heure en radian du coucher de soleil
148 c -omega est donc l'heure en radian de lever du soleil
149  real omegadeb, omegafin
150  real zfrac1, zfrac2, z1_mu, z2_mu
151  real lat_sun ! declinaison en radian
152  real lon_sun ! longitude solaire en radian
153  real latr ! latitude du pt de grille en radian
154 c================================================================
155 c
156  pi_local = 4.0 * atan(1.0)
157  deux_pi_local = 2.0 * pi_local
158  incl=r_incl * pi_local / 180.
159 c
160  lon_sun = longi * pi_local / 180.0
161  lat_sun = asin(sin(lon_sun)*sin(incl) )
162 c
163  gmtime1=gmtime*86400.
164  gmtime2=gmtime*86400.+pdtrad
165 c
166  DO i = 1, klon
167 c
168  latr = lat(i) * pi_local / 180.
169 c
170 c--pose probleme quand lat=+/-90 degres
171 c
172 c omega = -TAN(latr)*TAN(lat_sun)
173 c omega = ACOS(omega)
174 c IF (latr.GE.(pi_local/2.+lat_sun)
175 c . .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN
176 c omega = 0.0 ! nuit polaire
177 c ENDIF
178 c IF (latr.GE.(pi_local/2.-lat_sun)
179 c . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
180 c omega = pi_local ! journee polaire
181 c ENDIF
182 c
183 c--remplace par cela (le cas par defaut est different)
184 c
185  omega=0.0 !--nuit polaire
186  IF (latr.GE.(pi_local/2.-lat_sun)
187  . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
188  omega = pi_local ! journee polaire
189  ENDIF
190  IF (latr.LT.(pi_local/2.+lat_sun).AND.
191  . latr.GT.(-pi_local/2.+lat_sun).AND.
192  . latr.LT.(pi_local/2.-lat_sun).AND.
193  . latr.GT.(-pi_local/2.-lat_sun)) THEN
194  omega = -tan(latr)*tan(lat_sun)
195  omega = acos(omega)
196  ENDIF
197 c
198  omega1 = gmtime1 + long(i)*86400.0/360.0
199  omega1 = omega1 / 86400.0*deux_pi_local
200  omega1 = mod(omega1+deux_pi_local, deux_pi_local)
201  omega1 = omega1 - pi_local
202 c
203  omega2 = gmtime2 + long(i)*86400.0/360.0
204  omega2 = omega2 / 86400.0*deux_pi_local
205  omega2 = mod(omega2+deux_pi_local, deux_pi_local)
206  omega2 = omega2 - pi_local
207 c
208  IF (omega1.LE.omega2) THEN !--on est dans la meme journee locale
209 c
210  IF (omega2.LE.-omega .OR. omega1.GE.omega
211  . .OR. omega.LT.1e-5) THEN !--nuit
212  frac(i)=0.0
213  pmu0(i)=0.0
214  ELSE !--jour+nuit/jour
215  omegadeb=max(-omega,omega1)
216  omegafin=min(omega,omega2)
217  frac(i)=(omegafin-omegadeb)/(omega2-omega1)
218  pmu0(i)=sin(latr)*sin(lat_sun) +
219  . cos(latr)*cos(lat_sun)*
220  . (sin(omegafin)-sin(omegadeb))/
221  . (omegafin-omegadeb)
222  ENDIF
223 c
224  ELSE !---omega1 GT omega2 -- a cheval sur deux journees
225 c
226 c-------------------entre omega1 et pi
227  IF (omega1.GE.omega) THEN !--nuit
228  zfrac1=0.0
229  z1_mu =0.0
230  ELSE !--jour+nuit
231  omegadeb=max(-omega,omega1)
232  omegafin=omega
233  zfrac1=omegafin-omegadeb
234  z1_mu =sin(latr)*sin(lat_sun) +
235  . cos(latr)*cos(lat_sun)*
236  . (sin(omegafin)-sin(omegadeb))/
237  . (omegafin-omegadeb)
238  ENDIF
239 c---------------------entre -pi et omega2
240  IF (omega2.LE.-omega) THEN !--nuit
241  zfrac2=0.0
242  z2_mu =0.0
243  ELSE !--jour+nuit
244  omegadeb=-omega
245  omegafin=min(omega,omega2)
246  zfrac2=omegafin-omegadeb
247  z2_mu =sin(latr)*sin(lat_sun) +
248  . cos(latr)*cos(lat_sun)*
249  . (sin(omegafin)-sin(omegadeb))/
250  . (omegafin-omegadeb)
251 c
252  ENDIF
253 c-----------------------moyenne
254  frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1)
255  pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/max(zfrac1+zfrac2,1.e-10)
256 c
257  ENDIF !---comparaison omega1 et omega2
258 c
259  ENDDO
260 c
261  END
262 c===================================================================
263  SUBROUTINE zenith (longi, gmtime, lat, long,
264  s pmu0, fract)
265  USE dimphy
266  IMPLICIT none
267 c
268 c Auteur(s): Z.X. Li (LMD/ENS)
269 c
270 c Objet: calculer le cosinus de l'angle zenithal du soleil en
271 c connaissant la declinaison du soleil, la latitude et la
272 c longitude du point sur la terre, et le temps universel
273 c
274 c Arguments d'entree:
275 c longi : declinaison du soleil (en degres)
276 c gmtime : temps universel en second qui varie entre 0 et 86400
277 c lat : latitude en degres
278 c long : longitude en degres
279 c Arguments de sortie:
280 c pmu0 : cosinus de l'angle zenithal
281 c
282 c====================================================================
283 cym#include "dimensions.h"
284 cym#include "dimphy.h"
285 #include "YOMCST.h"
286 c====================================================================
287  REAL longi, gmtime
288  REAL lat(klon), long(klon), pmu0(klon), fract(klon)
289 c=====================================================================
290  INTEGER n
291  REAL zpi, zpir, omega, zgmtime
292  REAL incl, lat_sun, lon_sun
293 c----------------------------------------------------------------------
294  zpi = 4.0*atan(1.0)
295  zpir = zpi / 180.0
296  zgmtime=gmtime*86400.
297 c
298  incl=r_incl * zpir
299 c
300  lon_sun = longi * zpir
301  lat_sun = asin(sin(lon_sun)*sin(incl) )
302 c
303 c--initialisation a la nuit
304 c
305  DO n =1, klon
306  pmu0(n)=0.
307  fract(n)=0.0
308  ENDDO
309 c
310 c 1 degre en longitude = 240 secondes en temps
311 c
312  DO n = 1, klon
313  omega = zgmtime + long(n)*86400.0/360.0
314  omega = omega / 86400.0 * 2.0 * zpi
315  omega = mod(omega + 2.0 * zpi, 2.0 * zpi)
316  omega = omega - zpi
317  pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun)
318  . + cos(lat(n)*zpir) * cos(lat_sun)
319  . * cos(omega)
320  pmu0(n) = max(pmu0(n), 0.0)
321  IF (pmu0(n).GT.1.e-6) fract(n)=1.0
322  ENDDO
323 c
324  RETURN
325  END