/[lmdze]/trunk/phylmd/orbite.f
ViewVC logotype

Diff of /trunk/phylmd/orbite.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 20 by guez, Wed Oct 15 16:19:57 2008 UTC
# Line 121  contains Line 121  contains
121      use dimphy, only: klon      use dimphy, only: klon
122      use YOMCST, only: r_incl      use YOMCST, only: r_incl
123      use phyetat0_m, only: rlat, rlon      use phyetat0_m, only: rlat, rlon
124        use comconst, only: pi
125    
126      ! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et      ! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et
127      ! "angle" de Z.X. Li      ! "angle" de Z.X. Li
# Line 132  contains Line 133  contains
133      ! moyennes de "pmu0" et non des valeurs instantanées.      ! moyennes de "pmu0" et non des valeurs instantanées.
134      ! Du coup "frac" prend toutes les valeurs entre 0 et 1.      ! Du coup "frac" prend toutes les valeurs entre 0 et 1.
135    
136      ! Date : premiere version le 13 decembre 1994      ! Date : première version le 13 decembre 1994
137      !          revu pour  GCM  le 30 septembre 1996      !        revu pour  GCM  le 30 septembre 1996
138    
139      real, intent(in):: longi      real, intent(in):: longi
140      ! (longitude vraie de la terre dans son plan solaire a partir de      ! (longitude vraie de la terre dans son plan solaire a partir de
141      ! l'equinoxe de printemps (degre))      ! l'equinoxe de printemps) (in degrees)
142    
143      real, intent(in):: gmtime ! temps universel en fraction de jour      real, intent(in):: gmtime ! temps universel en fraction de jour
144      real, intent(in):: pdtrad ! pas de temps du rayonnement (secondes)      real, intent(in):: pdtrad ! pas de temps du rayonnement (secondes)
145    
146      real, intent(out):: pmu0(klon)      real, intent(out):: pmu0(klon)
147      ! (cosinus de l'angle zenithal moyen entre gmtime et gmtime+pdtrad)      ! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
148    
149      real, intent(out), optional:: frac(klon)      real, intent(out), optional:: frac(klon)
150      ! (ensoleillement moyen entre gmtime et gmtime+pdtrad)      ! (ensoleillement moyen entre gmtime et gmtime+pdtrad)
# Line 152  contains Line 153  contains
153    
154      integer i      integer i
155      real gmtime1, gmtime2      real gmtime1, gmtime2
156      real pi_local, deux_pi_local, incl      real deux_pi, incl
157    
158      real omega1, omega2, omega      real omega1, omega2, omega
159      ! omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi.      ! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.
160      ! omega : heure en radian du coucher de soleil      ! omega : heure en radians du coucher de soleil
161      ! -omega est donc l'heure en radian de lever du soleil      ! -omega est donc l'heure en radians de lever du soleil
162    
163      real omegadeb, omegafin      real omegadeb, omegafin
164      real zfrac1, zfrac2, z1_mu, z2_mu      real zfrac1, zfrac2, z1_mu, z2_mu
165      real lat_sun          ! declinaison en radian      real lat_sun          ! déclinaison en radians
166      real lon_sun          ! longitude solaire en radian      real lon_sun          ! longitude solaire en radians
167      real latr             ! latitude du pt de grille en radian      real latr             ! latitude du point de grille en radians
168    
169      !----------------------------------------------------------------------      !----------------------------------------------------------------------
170    
171      pi_local = 4.0 * ATAN(1.0)      deux_pi = 2 * pi
172      deux_pi_local = 2.0 * pi_local      incl=R_incl * pi / 180.
     incl=R_incl * pi_local / 180.  
173    
174      lon_sun = longi * pi_local / 180.0      lon_sun = longi * pi / 180.
175      lat_sun = ASIN (SIN(lon_sun)*SIN(incl) )      lat_sun = ASIN(SIN(lon_sun)*SIN(incl))
176    
177      gmtime1=gmtime*86400.      gmtime1=gmtime*86400.
178      gmtime2=gmtime*86400.+pdtrad      gmtime2=gmtime*86400.+pdtrad
179    
180      DO i = 1, klon      DO i = 1, klon
181         latr = rlat(i) * pi_local / 180.         latr = rlat(i) * pi / 180.
182         omega=0.0  !--nuit polaire         omega=0.0  !--nuit polaire
183         IF (latr >= (pi_local/2.-lat_sun) &         IF (latr >= (pi/2.-lat_sun) &
184              .OR. latr <= (-pi_local/2.-lat_sun)) THEN              .OR. latr <= (-pi/2.-lat_sun)) THEN
185            omega = pi_local  ! journee polaire            omega = pi  ! journee polaire
186         ENDIF         ENDIF
187         IF (latr < (pi_local/2.+lat_sun).AND. &         IF (latr < (pi/2.+lat_sun).AND. &
188              latr.GT.(-pi_local/2.+lat_sun).AND. &              latr.GT.(-pi/2.+lat_sun).AND. &
189              latr < (pi_local/2.-lat_sun).AND. &              latr < (pi/2.-lat_sun).AND. &
190              latr.GT.(-pi_local/2.-lat_sun)) THEN              latr.GT.(-pi/2.-lat_sun)) THEN
191            omega = -TAN(latr)*TAN(lat_sun)            omega = -TAN(latr)*TAN(lat_sun)
192            omega = ACOS(omega)            omega = ACOS(omega)
193         ENDIF         ENDIF
194    
195         omega1 = gmtime1 + rlon(i)*86400.0/360.0         omega1 = gmtime1 + rlon(i)*86400.0/360.0
196         omega1 = omega1 / 86400.0*deux_pi_local         omega1 = omega1 / 86400.0*deux_pi
197         omega1 = MOD (omega1+deux_pi_local, deux_pi_local)         omega1 = MOD (omega1+deux_pi, deux_pi)
198         omega1 = omega1 - pi_local         omega1 = omega1 - pi
199    
200         omega2 = gmtime2 + rlon(i)*86400.0/360.0         omega2 = gmtime2 + rlon(i)*86400.0/360.0
201         omega2 = omega2 / 86400.0*deux_pi_local         omega2 = omega2 / 86400.0*deux_pi
202         omega2 = MOD (omega2+deux_pi_local, deux_pi_local)         omega2 = MOD (omega2+deux_pi, deux_pi)
203         omega2 = omega2 - pi_local         omega2 = omega2 - pi
204    
205         test_omega12: IF (omega1 <= omega2) THEN           test_omega12: IF (omega1 <= omega2) THEN  
206            ! on est dans la meme journee locale            ! on est dans la meme journee locale
# Line 247  contains Line 249  contains
249            ENDIF            ENDIF
250            !-----------------------moyenne            !-----------------------moyenne
251            if (present(frac)) &            if (present(frac)) &
252                 frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1)                 frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi-omega1)
253            pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2, 1.E-10)            pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2, 1.E-10)
254         ENDIF test_omega12         ENDIF test_omega12
255      ENDDO      ENDDO
Legend:
Removed from v.3  
changed lines
  Added in v.20
  ViewVC Help
Powered by ViewVC 1.1.21