source: LMDZ5/trunk/libf/phylmd/orbite.F90 @ 2306

Last change on this file since 2306 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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