source: lmdz_wrf/trunk/WRFV3/lmdz/orbite.F90 @ 1419

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