source: LMDZ4/tags/pre_parallel/libf/phylmd/orbite.F @ 730

Last change on this file since 730 was 730, checked in by (none), 18 years ago

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