source: lmdz_wrf/trunk/WRFV3/lmdz/orbite_mod.F90 @ 1465

Last change on this file since 1465 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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