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

Diff of /trunk/phylmd/zenang.f

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

revision 117 by guez, Tue May 13 17:23:16 2014 UTC revision 118 by guez, Thu Dec 18 17:30:24 2014 UTC
# Line 4  module zenang_m Line 4  module zenang_m
4    
5  contains  contains
6    
7    SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)    SUBROUTINE zenang(longi, gmtime, pdtrad, mu0, frac)
8    
9      USE dimphy, ONLY : klon      USE dimphy, ONLY: klon
10      USE yomcst, ONLY : r_incl      USE yomcst, ONLY: r_incl
11      USE phyetat0_m, ONLY : rlat, rlon      USE phyetat0_m, ONLY: rlat, rlon
12      use nr_util, only: assert, pi      use nr_util, only: assert, pi
13    
14      ! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et      ! Author: O. Boucher (LMD/CNRS), d'après les routines "zenith" et
15      ! "angle" de Z.X. Li      ! "angle" de Z.X. Li
16    
17      ! Calcule les valeurs moyennes du cos de l'angle zénithal et      ! Date : première version le 13 décembre 1994, revu pour GCM le 30
18      ! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la      ! septembre 1996
     ! déclinaison, la latitude et la longitude.  
     ! Différent de la routine "angle" en ce sens que "zenang" fournit des  
     ! moyennes de "pmu0" et non des valeurs instantanées.  
     ! Du coup "frac" prend toutes les valeurs entre 0 et 1.  
   
     ! Date : première version le 13 decembre 1994  
     !        revu pour  GCM  le 30 septembre 1996  
19    
20      REAL, INTENT (IN):: longi      ! Calcule les valeurs moyennes du cos de l'angle zénithal et
21      ! (longitude vraie de la terre dans son plan solaire a partir de      ! l'ensoleillement moyen entre "gmtime" et "gmtime+pdtrad"
22      ! l'equinoxe de printemps) (in degrees)      ! connaissant la déclinaison, la latitude et la longitude.
23        ! Différent de la routine "angle" en ce sens que "zenang" fournit
24      REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour      ! des moyennes de "mu0" et non des valeurs instantanées.  Du coup
25      REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)      ! "frac" prend toutes les valeurs entre 0 et 1.
26    
27        REAL, INTENT(IN):: longi
28        ! longitude vraie de la terre dans son plan solaire à partir de
29        ! l'équinoxe de printemps (in degrees)
30    
31        REAL, INTENT(IN):: gmtime ! temps universel en fraction de jour
32        REAL, INTENT(IN):: pdtrad ! pas de temps du rayonnement (s)
33    
34      REAL, INTENT (OUT):: pmu0(:) ! (klon)      REAL, INTENT(OUT):: mu0(:) ! (klon)
35      ! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")      ! cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad"
36    
37      REAL, INTENT (OUT), OPTIONAL:: frac(:) ! (klon)      REAL, INTENT(OUT), OPTIONAL:: frac(:) ! (klon)
38      ! (ensoleillement moyen entre gmtime et gmtime+pdtrad)      ! ensoleillement moyen entre gmtime et gmtime+pdtrad
39    
40      ! Variables local to the procedure:      ! Local:
41    
42      INTEGER i      INTEGER i
43      REAL gmtime1, gmtime2      REAL gmtime1, gmtime2
44      REAL deux_pi      REAL deux_pi
45    
46      REAL omega1, omega2, omega      REAL omega1, omega2 ! temps 1 et 2 exprimés en radians avec 0 à midi
47      ! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.  
48      ! omega : heure en radians du coucher de soleil      REAL omega ! heure en rad du coucher de soleil
49      ! -omega est donc l'heure en radians de lever du soleil      ! - omega est donc l'heure en rad de lever du soleil
50    
51      REAL omegadeb, omegafin      REAL omegadeb, omegafin
52      REAL zfrac1, zfrac2, z1_mu, z2_mu      REAL zfrac1, zfrac2, z1_mu, z2_mu
# Line 55  contains Line 55  contains
55    
56      !----------------------------------------------------------------------      !----------------------------------------------------------------------
57    
58      if (present(frac)) call assert((/size(pmu0), size(frac)/) == klon, &      if (present(frac)) call assert((/size(mu0), size(frac)/) == klon, "zenang")
          "zenang")  
59            
60      deux_pi = 2*pi      deux_pi = 2*pi
61    
62      lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))      lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))
63      ! Capderou (2003 #784, équation 4.49)      ! Capderou (2003 784, equation 4.49)
64    
65      gmtime1 = gmtime*86400.      gmtime1 = gmtime*86400.
66      gmtime2 = gmtime*86400. + pdtrad      gmtime2 = gmtime*86400. + pdtrad
67    
68      DO i = 1, klon      DO i = 1, klon
69         latr = rlat(i)*pi/180.         latr = rlat(i)*pi/180.
70         omega = 0.0 !--nuit polaire         omega = 0.0 ! nuit polaire
71         IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN         IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN
72            omega = pi ! journee polaire            omega = pi ! journée polaire
73         END IF         END IF
74         IF (latr<(pi/2.+lat_sun) .AND. latr>(-pi/2.+lat_sun) .AND. &         IF (latr<(pi/2.+lat_sun) .AND. latr>(-pi/2.+lat_sun) .AND. &
75              latr<(pi/2.-lat_sun) .AND. latr>(-pi/2.-lat_sun)) THEN              latr<(pi/2.-lat_sun) .AND. latr>(-pi/2.-lat_sun)) THEN
# Line 88  contains Line 87  contains
87         omega2 = mod(omega2+deux_pi, deux_pi)         omega2 = mod(omega2+deux_pi, deux_pi)
88         omega2 = omega2 - pi         omega2 = omega2 - pi
89    
90         TEST_OMEGA12: IF (omega1<=omega2) THEN         IF (omega1<=omega2) THEN
91            ! on est dans la meme journee locale            ! on est dans la meme journee locale
92            IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN            IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN
93               ! nuit               ! nuit
94               IF (present(frac)) frac(i) = 0.0               IF (present(frac)) frac(i) = 0.0
95               pmu0(i) = 0.0               mu0(i) = 0.0
96            ELSE            ELSE
97               ! jour + nuit / jour               ! jour + nuit / jour
98               omegadeb = max(-omega, omega1)               omegadeb = max(-omega, omega1)
99               omegafin = min(omega, omega2)               omegafin = min(omega, omega2)
100               IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)               IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)
101               pmu0(i) = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &               mu0(i) = sin(latr) * sin(lat_sun) + cos(latr) * cos(lat_sun) &
102                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)                    * (sin(omegafin) - sin(omegadeb)) / (omegafin - omegadeb)
103            END IF            END IF
104         ELSE TEST_OMEGA12         ELSE
105            !---omega1 GT omega2 -- a cheval sur deux journees            ! omega1 > omega2, à cheval sur deux journées
106            !-------------------entre omega1 et pi            ! entre omega1 et pi
107            IF (omega1>=omega) THEN !--nuit            IF (omega1>=omega) THEN ! nuit
108               zfrac1 = 0.0               zfrac1 = 0.0
109               z1_mu = 0.0               z1_mu = 0.0
110            ELSE !--jour+nuit            ELSE ! jour+nuit
111               omegadeb = max(-omega, omega1)               omegadeb = max(-omega, omega1)
112               omegafin = omega               omegafin = omega
113               zfrac1 = omegafin - omegadeb               zfrac1 = omegafin - omegadeb
114               z1_mu = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &               z1_mu = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &
115                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)
116            END IF            END IF
117            !---------------------entre -pi et omega2            ! entre -pi et omega2
118            IF (omega2<=-omega) THEN !--nuit            IF (omega2<=-omega) THEN ! nuit
119               zfrac2 = 0.0               zfrac2 = 0.0
120               z2_mu = 0.0               z2_mu = 0.0
121            ELSE !--jour+nuit            ELSE ! jour+nuit
122               omegadeb = -omega               omegadeb = -omega
123               omegafin = min(omega, omega2)               omegafin = min(omega, omega2)
124               zfrac2 = omegafin - omegadeb               zfrac2 = omegafin - omegadeb
# Line 127  contains Line 126  contains
126                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)                    omegafin)-sin(omegadeb))/(omegafin-omegadeb)
127    
128            END IF            END IF
129            !-----------------------moyenne            ! moyenne
130            IF (present(frac)) frac(i) = (zfrac1+zfrac2)/ &            IF (present(frac)) frac(i) = (zfrac1+zfrac2)/ (omega2+deux_pi-omega1)
131                 (omega2+deux_pi-omega1)            mu0(i) = (zfrac1*z1_mu+zfrac2*z2_mu)/max(zfrac1+zfrac2, 1.E-10)
132            pmu0(i) = (zfrac1*z1_mu+zfrac2*z2_mu)/max(zfrac1+zfrac2, 1.E-10)         END IF
        END IF TEST_OMEGA12  
133      END DO      END DO
134    
135    END SUBROUTINE zenang    END SUBROUTINE zenang
Legend:
Removed from v.117  
changed lines
  Added in v.118
  ViewVC Help
Powered by ViewVC 1.1.21