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

Diff of /trunk/Sources/phylmd/orbite.f

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

revision 20 by guez, Wed Oct 15 16:19:57 2008 UTC revision 22 by guez, Fri Jul 31 15:18:47 2009 UTC
# Line 1  Line 1 
1  module orbite_m  MODULE orbite_m
2    
3    ! From phylmd/orbite.F, v 1.1.1.1 2004/05/19 12:53:08    ! From phylmd/orbite.F, v 1.1.1.1 2004/05/19 12:53:08
4    
5    IMPLICIT none    IMPLICIT NONE
6    
7  contains  CONTAINS
8    
9    SUBROUTINE orbite(xjour, longi, dist)    SUBROUTINE orbite(xjour, longi, dist)
10    
11      use YOMCST      USE yomcst, ONLY : r_ecc, r_peri
12        use numer_rec, only: pi
13    
14      ! Auteur(s): Z.X. Li (LMD/CNRS) Date: 1993/08/18 Pour un jour      ! Auteur(s): Z.X. Li (LMD/CNRS)
15      ! donné, calcule la longitude vraie de la Terre (par rapport au      ! Date: 1993/08/18
16      ! point vernal, 21 mars) dans son orbite solaire. Calcule aussi      ! Pour un jour donné, calcule la longitude vraie de la Terre (par
17        ! rapport au point vernal, 21 mars) dans son orbite solaire. Calcule aussi
18      ! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.      ! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.
19    
20      REAL, intent(in):: xjour ! jour de l'annee a compter du 1er janvier      REAL, INTENT (IN):: xjour ! jour de l'année à compter du premier janvier
21    
22      real, intent(out):: longi      REAL, INTENT (OUT):: longi
23      ! (longitude vraie de la Terre dans son orbite solaire, par      ! longitude vraie de la Terre dans son orbite solaire, par
24      ! rapport au point vernal (21 mars), en degrés)      ! rapport au point vernal (21 mars), en degrés
25    
26      real, intent(out), optional:: dist      REAL, INTENT (OUT), OPTIONAL:: dist
27      ! (distance terre-soleil (par rapport a la moyenne))      ! distance terre-soleil (par rapport a la moyenne)
28    
29      ! Variables locales      ! Variables locales
30      REAL pir, xl, xllp, xee, xse, xlam, dlamm, anm, ranm, anv, ranv      REAL pir, xl, xllp, xee, xse, xlam, anm, ranm, ranv
31    
32      !----------------------------------------------------------------------      !----------------------------------------------------------------------
33    
34      pir = 4.0*ATAN(1.0) / 180.0      pir = pi / 180.
35      xl=R_peri+180.0      xl = r_peri + 180.
36      xllp=xl*pir      xllp = xl * pir
37      xee=R_ecc*R_ecc      xee = r_ecc * r_ecc
38      xse=SQRT(1.0-xee)      xse = sqrt(1. - xee)
39      xlam = (R_ecc/2.0+R_ecc*xee/8.0)*(1.0+xse)*SIN(xllp) &      xlam = (r_ecc / 2 + r_ecc * xee / 8.) * (1. + xse) * sin(xllp) &
40           - xee/4.0*(0.5+xse)*SIN(2.0*xllp) &           - xee / 4. * (0.5 + xse) * sin(2.*xllp) &
41           + R_ecc*xee/8.0*(1.0/3.0+xse)*SIN(3.0*xllp)           + r_ecc * xee / 8. * (1. / 3. + xse) * sin(3.*xllp)
42      xlam=2.0*xlam/pir      xlam = 2. * xlam / pir
43      dlamm=xlam+(xjour-81.0)      anm = xlam + (xjour - 81.) - xl
44      anm=dlamm-xl      ranm = anm * pir
45      ranm=anm*pir      xee = xee * r_ecc
46      xee=xee*R_ecc      ranv = ranm + (2. * r_ecc - xee / 4.) * sin(ranm) + &
47      ranv=ranm+(2.0*R_ecc-xee/4.0)*SIN(ranm) &           5. / 4. * r_ecc * r_ecc * sin(2 * ranm) &
48           +5.0/4.0*R_ecc*R_ecc*SIN(2.0*ranm) &           + 13. / 12. * xee * sin(3.*ranm)
49           +13.0/12.0*xee*SIN(3.0*ranm)  
50        longi = ranv / pir + xl
51      anv=ranv/pir  
52      longi=anv+xl      IF (present(dist)) then
53           dist = (1 - r_ecc*r_ecc) &
54      if (present(dist)) dist &              / (1 + r_ecc*cos(pir*(longi - (r_peri + 180.))))
55           = (1-R_ecc*R_ecc) /(1+R_ecc*COS(pir*(longi-(R_peri+180.0))))      end IF
56    
57    END SUBROUTINE orbite    END SUBROUTINE orbite
58    
59    !*****************************************************************    !*****************************************************************
60    
   SUBROUTINE angle(longi, lati, frac, muzero)  
   
     use dimphy, only: klon  
     use YOMCST, only: R_incl  
   
     ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18  
     ! Calcule la durée d'ensoleillement pour un jour et la hauteur du  
     ! soleil (cosinus de l'angle zénithal) moyennée sur la journee.  
   
     ! Arguments:  
     ! longi----INPUT-R- la longitude vraie de la terre dans son plan  
     !                   solaire a partir de l'equinoxe de printemps (degre)  
     ! lati-----INPUT-R- la latitude d'un point sur la terre (degre)  
     ! frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee  
     !                   par 24 heures (unite en fraction de 0 a 1)  
     ! muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur  
     !                   la journee (0 a 1)  
   
     REAL longi  
     REAL lati(klon), frac(klon), muzero(klon)  
     REAL lat, omega, lon_sun, lat_sun  
     REAL pi_local, incl  
     INTEGER i  
   
     !----------------------------------------------------------------------  
   
     pi_local = 4.0 * ATAN(1.0)  
     incl=R_incl * pi_local / 180.  
   
     lon_sun = longi * pi_local / 180.0  
     lat_sun = ASIN (sin(lon_sun)*SIN(incl) )  
   
     DO i = 1, klon  
        lat = lati(i) * pi_local / 180.0  
   
        IF ( lat .GE. (pi_local/2.+lat_sun) &  
             .OR. lat <= (-pi_local/2.+lat_sun)) THEN  
           omega = 0.0   ! nuit polaire  
        ELSE IF ( lat.GE.(pi_local/2.-lat_sun) &  
             .OR. lat <= (-pi_local/2.-lat_sun)) THEN  
           omega = pi_local   ! journee polaire  
        ELSE  
           omega = -TAN(lat)*TAN(lat_sun)  
           omega = ACOS (omega)  
        ENDIF  
   
        frac(i) = omega / pi_local  
   
        IF (omega .GT. 0.0) THEN  
           muzero(i) = SIN(lat)*SIN(lat_sun) &  
                + COS(lat)*COS(lat_sun)*SIN(omega) / omega  
        ELSE  
           muzero(i) = 0.0  
        ENDIF  
     ENDDO  
   
   END SUBROUTINE angle  
   
   !*****************************************************************  
   
61    SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)    SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)
62    
63      use dimphy, only: klon      USE dimphy, ONLY : klon
64      use YOMCST, only: r_incl      USE yomcst, ONLY : r_incl
65      use phyetat0_m, only: rlat, rlon      USE phyetat0_m, ONLY : rlat, rlon
66      use comconst, only: pi      USE comconst, ONLY : pi
67        use numer_rec, only: assert
68    
69      ! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et      ! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et
70      ! "angle" de Z.X. Li      ! "angle" de Z.X. Li
# Line 136  contains Line 79  contains
79      ! Date : première version le 13 decembre 1994      ! Date : première version le 13 decembre 1994
80      !        revu pour  GCM  le 30 septembre 1996      !        revu pour  GCM  le 30 septembre 1996
81    
82      real, intent(in):: longi      REAL, INTENT (IN):: longi
83      ! (longitude vraie de la terre dans son plan solaire a partir de      ! (longitude vraie de la terre dans son plan solaire a partir de
84      ! l'equinoxe de printemps) (in degrees)      ! l'equinoxe de printemps) (in degrees)
85    
86      real, intent(in):: gmtime ! temps universel en fraction de jour      REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour
87      real, intent(in):: pdtrad ! pas de temps du rayonnement (secondes)      REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)
88    
89      real, intent(out):: pmu0(klon)      REAL, INTENT (OUT):: pmu0(:) ! (klon)
90      ! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")      ! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")
91    
92      real, intent(out), optional:: frac(klon)      REAL, INTENT (OUT), OPTIONAL:: frac(:) ! (klon)
93      ! (ensoleillement moyen entre gmtime et gmtime+pdtrad)      ! (ensoleillement moyen entre gmtime et gmtime+pdtrad)
94    
95      ! Variables local to the procedure:      ! Variables local to the procedure:
96    
97      integer i      INTEGER i
98      real gmtime1, gmtime2      REAL gmtime1, gmtime2
99      real deux_pi, incl      REAL deux_pi
100    
101      real omega1, omega2, omega      REAL omega1, omega2, omega
102      ! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.      ! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.
103      ! omega : heure en radians du coucher de soleil      ! omega : heure en radians du coucher de soleil
104      ! -omega est donc l'heure en radians de lever du soleil      ! -omega est donc l'heure en radians de lever du soleil
105    
106      real omegadeb, omegafin      REAL omegadeb, omegafin
107      real zfrac1, zfrac2, z1_mu, z2_mu      REAL zfrac1, zfrac2, z1_mu, z2_mu
108      real lat_sun          ! déclinaison en radians      REAL lat_sun ! déclinaison en radians
109      real lon_sun          ! longitude solaire en radians      REAL latr ! latitude du point de grille en radians
     real latr             ! latitude du point de grille en radians  
110    
111      !----------------------------------------------------------------------      !----------------------------------------------------------------------
112    
113      deux_pi = 2 * pi      call assert((/size(pmu0), size(frac)/) == klon, "zenang")
114      incl=R_incl * pi / 180.  
115        deux_pi = 2*pi
116    
117      lon_sun = longi * pi / 180.      lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
118      lat_sun = ASIN(SIN(lon_sun)*SIN(incl))      ! Capderou (2003 #784, équation 4.49)
119    
120      gmtime1=gmtime*86400.      gmtime1 = gmtime*86400.
121      gmtime2=gmtime*86400.+pdtrad      gmtime2 = gmtime*86400. + pdtrad
122    
123      DO i = 1, klon      DO i = 1, klon
124         latr = rlat(i) * pi / 180.         latr = rlat(i)*pi/180.
125         omega=0.0  !--nuit polaire         omega = 0.0 !--nuit polaire
126         IF (latr >= (pi/2.-lat_sun) &         IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN
127              .OR. latr <= (-pi/2.-lat_sun)) THEN            omega = pi ! journee polaire
128            omega = pi  ! journee polaire         END IF
129         ENDIF         IF (latr<(pi/2.+lat_sun) .AND. latr>(-pi/2.+lat_sun) .AND. &
130         IF (latr < (pi/2.+lat_sun).AND. &              latr<(pi/2.-lat_sun) .AND. latr>(-pi/2.-lat_sun)) THEN
131              latr.GT.(-pi/2.+lat_sun).AND. &            omega = -tan(latr)*tan(lat_sun)
132              latr < (pi/2.-lat_sun).AND. &            omega = acos(omega)
133              latr.GT.(-pi/2.-lat_sun)) THEN         END IF
           omega = -TAN(latr)*TAN(lat_sun)  
           omega = ACOS(omega)  
        ENDIF  
134    
135         omega1 = gmtime1 + rlon(i)*86400.0/360.0         omega1 = gmtime1 + rlon(i)*86400.0/360.0
136         omega1 = omega1 / 86400.0*deux_pi         omega1 = omega1/86400.0*deux_pi
137         omega1 = MOD (omega1+deux_pi, deux_pi)         omega1 = mod(omega1+deux_pi, deux_pi)
138         omega1 = omega1 - pi         omega1 = omega1 - pi
139    
140         omega2 = gmtime2 + rlon(i)*86400.0/360.0         omega2 = gmtime2 + rlon(i)*86400.0/360.0
141         omega2 = omega2 / 86400.0*deux_pi         omega2 = omega2/86400.0*deux_pi
142         omega2 = MOD (omega2+deux_pi, deux_pi)         omega2 = mod(omega2+deux_pi, deux_pi)
143         omega2 = omega2 - pi         omega2 = omega2 - pi
144    
145         test_omega12: IF (omega1 <= omega2) THEN           TEST_OMEGA12: IF (omega1<=omega2) THEN
146            ! on est dans la meme journee locale            ! on est dans la meme journee locale
147            IF (omega2 <= -omega .OR. omega1 >= omega .OR. omega < 1e-5) THEN            IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
148               ! nuit               ! nuit
149               if (present(frac)) frac(i)=0.0               IF (present(frac)) frac(i) = 0.0
150               pmu0(i)=0.0               pmu0(i) = 0.0
151            ELSE            ELSE
152               ! jour + nuit / jour               ! jour + nuit / jour
153               omegadeb=MAX(-omega, omega1)               omegadeb = max(-omega, omega1)
154               omegafin=MIN(omega, omega2)               omegafin = min(omega, omega2)
155               if (present(frac)) frac(i)=(omegafin-omegadeb)/(omega2-omega1)               IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)
156               pmu0(i)=SIN(latr)*SIN(lat_sun) +  &               pmu0(i) = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &
157                    COS(latr)*COS(lat_sun)* &                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)
158                    (SIN(omegafin)-SIN(omegadeb))/ &            END IF
159                    (omegafin-omegadeb)                 ELSE TEST_OMEGA12
           ENDIF  
        ELSE test_omega12  
160            !---omega1 GT omega2 -- a cheval sur deux journees            !---omega1 GT omega2 -- a cheval sur deux journees
161            !-------------------entre omega1 et pi            !-------------------entre omega1 et pi
162            IF (omega1 >= omega) THEN  !--nuit            IF (omega1>=omega) THEN !--nuit
163               zfrac1=0.0               zfrac1 = 0.0
164               z1_mu =0.0               z1_mu = 0.0
165            ELSE                       !--jour+nuit            ELSE !--jour+nuit
166               omegadeb=MAX(-omega, omega1)               omegadeb = max(-omega, omega1)
167               omegafin=omega               omegafin = omega
168               zfrac1=omegafin-omegadeb               zfrac1 = omegafin - omegadeb
169               z1_mu =SIN(latr)*SIN(lat_sun) + &               z1_mu = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &
170                    COS(latr)*COS(lat_sun)* &                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)
171                    (SIN(omegafin)-SIN(omegadeb))/ &            END IF
                   (omegafin-omegadeb)  
           ENDIF  
172            !---------------------entre -pi et omega2            !---------------------entre -pi et omega2
173            IF (omega2 <= -omega) THEN   !--nuit            IF (omega2<=-omega) THEN !--nuit
174               zfrac2=0.0               zfrac2 = 0.0
175               z2_mu =0.0               z2_mu = 0.0
176            ELSE                         !--jour+nuit            ELSE !--jour+nuit
177               omegadeb=-omega               omegadeb = -omega
178               omegafin=MIN(omega, omega2)               omegafin = min(omega, omega2)
179               zfrac2=omegafin-omegadeb               zfrac2 = omegafin - omegadeb
180               z2_mu =SIN(latr)*SIN(lat_sun) + &               z2_mu = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &
181                    COS(latr)*COS(lat_sun)* &                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)
                   (SIN(omegafin)-SIN(omegadeb))/ &  
                   (omegafin-omegadeb)  
182    
183            ENDIF            END IF
184            !-----------------------moyenne            !-----------------------moyenne
185            if (present(frac)) &            IF (present(frac)) frac(i) = (zfrac1+zfrac2)/ &
186                 frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi-omega1)                 (omega2+deux_pi-omega1)
187            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)
188         ENDIF test_omega12         END IF TEST_OMEGA12
189      ENDDO      END DO
190    
191    END SUBROUTINE zenang    END SUBROUTINE zenang
192    
193    !*****************************************************************  END MODULE orbite_m
   
   SUBROUTINE zenith (longi, gmtime, lat, long, pmu0, fract)  
   
     use dimphy, only: klon  
     use YOMCST, only: R_incl  
     use comconst, only: pi  
   
     ! Author: Z.X. Li (LMD/ENS)  
   
     ! Calcule le cosinus de l'angle zénithal du Soleil en connaissant  
     ! la déclinaison du Soleil, la latitude et la longitude du point  
     ! sur la Terre, et le temps universel.  
   
     REAL, intent(in):: longi ! déclinaison du soleil (en degrés)  
   
     REAL, intent(in):: gmtime  
     ! (temps universel en fraction de jour, entre 0 et 1)  
   
     REAL, intent(in):: lat(klon), long(klon) ! latitude et longitude en degres  
     REAL, intent(out):: pmu0(klon) ! cosinus de l'angle zenithal  
     REAL, intent(out):: fract(klon)  
   
     ! Variables local to the procedure:  
   
     INTEGER n  
     REAL zpir, omega  
     REAL lat_sun  
   
     !----------------------------------------------------------------------  
   
     zpir = pi / 180.0  
     lat_sun = ASIN (SIN(longi * zpir)*SIN(R_incl * zpir) )  
   
     ! Initialisation à la nuit:  
     pmu0 = 0.  
     fract = 0.  
   
     ! 1 degree in longitude = 240 s in time  
   
     DO n = 1, klon  
        omega = (gmtime + long(n) / 360.) * 2 * pi  
        omega = MOD(omega + 2 * pi, 2 * pi)  
        omega = omega - pi  
        pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun) &  
             + cos(lat(n)*zpir) * cos(lat_sun) &  
             * cos(omega)  
        pmu0(n) = MAX (pmu0(n), 0.)  
        IF (pmu0(n) > 1.E-6) fract(n)=1.0  
     ENDDO  
   
   END SUBROUTINE zenith  
   
 end module orbite_m  
Legend:
Removed from v.20  
changed lines
  Added in v.22
  ViewVC Help
Powered by ViewVC 1.1.21