source: LMDZ.3.3/trunk/libf/phylmd/orbite.F @ 5306

Last change on this file since 5306 was 2, checked in by lmdz, 25 years ago

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