source: LMDZ4/trunk/libf/phylmd/orbite.F @ 4077

Last change on this file since 4077 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
  • 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      USE dimphy
52      IMPLICIT none
53c======================================================================
54c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
55c Objet: Calculer la duree d'ensoleillement pour un jour et la hauteur
56c        du soleil (cosinus de l'angle zinithal) moyenne sur la journee
57c======================================================================
58c Arguments:
59c longi----INPUT-R- la longitude vraie de la terre dans son plan
60c                   solaire a partir de l'equinoxe de printemps (degre)
61c lati-----INPUT-R- la latitude d'un point sur la terre (degre)
62c frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee
63c                   par 24 heures (unite en fraction de 0 a 1)
64c muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur
65c                   la journee (0 a 1)
66c======================================================================
67cym#include "dimensions.h"
68cym#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
75c
76      pi_local = 4.0 * ATAN(1.0)
77      incl=R_incl * pi_local / 180.
78c
79      lon_sun = longi * pi_local / 180.0
80      lat_sun = ASIN (sin(lon_sun)*SIN(incl) )
81c
82      DO i = 1, klon
83      lat = lati(i) * pi_local / 180.0
84c
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
95c
96      frac(i) = omega / pi_local
97c
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
105c
106      RETURN
107      END
108c====================================================================
109      SUBROUTINE zenang(longi,gmtime,pdtrad,lat,long,
110     s                  pmu0,frac)
111      USE dimphy
112      IMPLICIT none
113c=============================================================
114c Auteur : O. Boucher (LMD/CNRS)
115c          d'apres les routines zenith et angle de Z.X. Li
116c Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
117c          et l'ensoleillement moyen entre gmtime1 et gmtime2
118c          connaissant la declinaison, la latitude et la longitude.
119c Rque   : Different de la routine angle en ce sens que zenang
120c          fournit des moyennes de pmu0 et non des valeurs
121c          instantanees, du coup frac prend toutes les valeurs
122c          entre 0 et 1.
123c Date   : premiere version le 13 decembre 1994
124c          revu pour  GCM  le 30 septembre 1996
125c===============================================================
126c longi : la longitude vraie de la terre dans son plan
127c                  solaire a partir de l'equinoxe de printemps (degre)
128c gmtime : temps universel en fraction de jour
129c pdtrad : pas de temps du rayonnement (secondes)
130c lat------INPUT : latitude en degres
131c long-----INPUT : longitude en degres
132c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
133c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
134c================================================================
135cym#include "dimensions.h"
136cym#include "dimphy.h"
137#include "YOMCST.h"
138c================================================================
139      real, intent(in):: longi, gmtime, pdtrad
140      real lat(klon), long(klon), pmu0(klon), frac(klon)
141c================================================================
142      integer i
143      real gmtime1, gmtime2
144      real pi_local, deux_pi_local, incl
145      real omega1, omega2, omega
146c omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.
147c omega : heure en radian du coucher de soleil
148c -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
154c================================================================
155c
156      pi_local = 4.0 * ATAN(1.0)
157      deux_pi_local = 2.0 * pi_local
158      incl=R_incl * pi_local / 180.
159c
160      lon_sun = longi * pi_local / 180.0
161      lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )
162c
163      gmtime1=gmtime*86400.
164      gmtime2=gmtime*86400.+pdtrad
165c
166      DO i = 1, klon
167c
168      latr = lat(i) * pi_local / 180.
169c
170c--pose probleme quand lat=+/-90 degres
171c
172c      omega = -TAN(latr)*TAN(lat_sun)
173c      omega = ACOS(omega)
174c      IF (latr.GE.(pi_local/2.+lat_sun)
175c     .    .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN
176c         omega = 0.0       ! nuit polaire
177c      ENDIF
178c      IF (latr.GE.(pi_local/2.-lat_sun)
179c     .          .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN
180c         omega = pi_local  ! journee polaire
181c      ENDIF
182c
183c--remplace par cela (le cas par defaut est different)
184c
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
197c
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
202c
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
207c
208      IF (omega1.LE.omega2) THEN  !--on est dans la meme journee locale
209c
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
223c
224      ELSE  !---omega1 GT omega2 -- a cheval sur deux journees
225c
226c-------------------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
239c---------------------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)
251c
252      ENDIF
253c-----------------------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)
256c
257      ENDIF   !---comparaison omega1 et omega2
258c
259      ENDDO
260c
261      END
262c===================================================================
263      SUBROUTINE zenith (longi, gmtime, lat, long,
264     s                   pmu0, fract)
265      USE dimphy
266      IMPLICIT none
267c
268c Auteur(s): Z.X. Li (LMD/ENS)
269c
270c Objet: calculer le cosinus de l'angle zenithal du soleil en
271c        connaissant la declinaison du soleil, la latitude et la
272c        longitude du point sur la terre, et le temps universel
273c
274c Arguments d'entree:
275c     longi  : declinaison du soleil (en degres)
276c     gmtime : temps universel en second qui varie entre 0 et 86400
277c     lat    : latitude en degres
278c     long   : longitude en degres
279c Arguments de sortie:
280c     pmu0   : cosinus de l'angle zenithal
281c
282c====================================================================
283cym#include "dimensions.h"
284cym#include "dimphy.h"
285#include "YOMCST.h"
286c====================================================================
287      REAL longi, gmtime
288      REAL lat(klon), long(klon), pmu0(klon), fract(klon)
289c=====================================================================
290      INTEGER n
291      REAL zpi, zpir, omega, zgmtime
292      REAL incl, lat_sun, lon_sun
293c----------------------------------------------------------------------
294      zpi = 4.0*ATAN(1.0)
295      zpir = zpi / 180.0
296      zgmtime=gmtime*86400.
297c
298      incl=R_incl * zpir
299c
300      lon_sun = longi * zpir
301      lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )
302c
303c--initialisation a la nuit
304c
305      DO n =1, klon
306        pmu0(n)=0.
307        fract(n)=0.0
308      ENDDO
309c
310c 1 degre en longitude = 240 secondes en temps
311c
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
323c
324      RETURN
325      END
Note: See TracBrowser for help on using the repository browser.