/[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

trunk/libf/phylmd/orbite.f90 revision 36 by guez, Thu Dec 2 17:11:04 2010 UTC trunk/phylmd/orbite.f revision 118 by guez, Thu Dec 18 17:30:24 2014 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, version 1.1.1.1 2004/05/19 12:53:08
4    
5    IMPLICIT NONE    IMPLICIT NONE
6    
# Line 8  CONTAINS Line 8  CONTAINS
8    
9    SUBROUTINE orbite(xjour, longi, dist)    SUBROUTINE orbite(xjour, longi, dist)
10    
     USE yomcst, ONLY : r_ecc, r_peri  
11      use nr_util, only: pi      use nr_util, only: pi
12        USE yomcst, ONLY : r_ecc, r_peri
13    
14      ! Auteur(s): Z.X. Li (LMD/CNRS)      ! Auteur(s): Z.X. Li (LMD/CNRS)
15      ! Date: 1993/08/18      ! Date: 1993/08/18
16    
17      ! Pour un jour donné, calcule la longitude vraie de la Terre (par      ! Pour un jour donné, calcule la longitude vraie de la Terre (par
18      ! rapport au point vernal, 21 mars) dans son orbite solaire. Calcule aussi      ! rapport au point vernal, 21 mars) dans son orbite
19      ! la distance Terre-Soleil, c'est-à-dire l'unité astronomique.      ! solaire. Calcule aussi la distance Terre-Soleil, c'est-à-dire
20        ! l'unité astronomique.
21    
22      REAL, INTENT (IN):: xjour ! jour de l'année à compter du premier janvier      REAL, INTENT (IN):: xjour ! jour de l'année à compter du premier janvier
23    
# Line 56  CONTAINS Line 58  CONTAINS
58    
59    END SUBROUTINE orbite    END SUBROUTINE orbite
60    
   !*****************************************************************  
   
   SUBROUTINE zenang(longi, gmtime, pdtrad, pmu0, frac)  
   
     USE dimphy, ONLY : klon  
     USE yomcst, ONLY : r_incl  
     USE phyetat0_m, ONLY : rlat, rlon  
     USE comconst, ONLY : pi  
     use nr_util, only: assert  
   
     ! Author : O. Boucher (LMD/CNRS), d'après les routines "zenith" et  
     ! "angle" de Z.X. Li  
   
     ! Calcule les valeurs moyennes du cos de l'angle zénithal et  
     ! l'ensoleillement moyen entre "gmtime1" et "gmtime2" connaissant la  
     ! 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  
   
     REAL, INTENT (IN):: longi  
     ! (longitude vraie de la terre dans son plan solaire a partir de  
     ! l'equinoxe de printemps) (in degrees)  
   
     REAL, INTENT (IN):: gmtime ! temps universel en fraction de jour  
     REAL, INTENT (IN):: pdtrad ! pas de temps du rayonnement (secondes)  
   
     REAL, INTENT (OUT):: pmu0(:) ! (klon)  
     ! (cosine of mean zenith angle between "gmtime" and "gmtime+pdtrad")  
   
     REAL, INTENT (OUT), OPTIONAL:: frac(:) ! (klon)  
     ! (ensoleillement moyen entre gmtime et gmtime+pdtrad)  
   
     ! Variables local to the procedure:  
   
     INTEGER i  
     REAL gmtime1, gmtime2  
     REAL deux_pi  
   
     REAL omega1, omega2, omega  
     ! omega1, omega2 : temps 1 et 2 exprimés en radians avec 0 à midi.  
     ! omega : heure en radians du coucher de soleil  
     ! -omega est donc l'heure en radians de lever du soleil  
   
     REAL omegadeb, omegafin  
     REAL zfrac1, zfrac2, z1_mu, z2_mu  
     REAL lat_sun ! déclinaison en radians  
     REAL latr ! latitude du point de grille en radians  
   
     !----------------------------------------------------------------------  
   
     if (present(frac)) call assert((/size(pmu0), size(frac)/) == klon, &  
          "zenang")  
       
     deux_pi = 2*pi  
   
     lat_sun = asin(sin(longi * pi / 180.) * sin(r_incl * pi / 180.))  
     ! Capderou (2003 #784, équation 4.49)  
   
     gmtime1 = gmtime*86400.  
     gmtime2 = gmtime*86400. + pdtrad  
   
     DO i = 1, klon  
        latr = rlat(i)*pi/180.  
        omega = 0.0 !--nuit polaire  
        IF (latr>=(pi/2.-lat_sun) .OR. latr<=(-pi/2.-lat_sun)) THEN  
           omega = pi ! journee polaire  
        END IF  
        IF (latr<(pi/2.+lat_sun) .AND. latr>(-pi/2.+lat_sun) .AND. &  
             latr<(pi/2.-lat_sun) .AND. latr>(-pi/2.-lat_sun)) THEN  
           omega = -tan(latr)*tan(lat_sun)  
           omega = acos(omega)  
        END IF  
   
        omega1 = gmtime1 + rlon(i)*86400.0/360.0  
        omega1 = omega1/86400.0*deux_pi  
        omega1 = mod(omega1+deux_pi, deux_pi)  
        omega1 = omega1 - pi  
   
        omega2 = gmtime2 + rlon(i)*86400.0/360.0  
        omega2 = omega2/86400.0*deux_pi  
        omega2 = mod(omega2+deux_pi, deux_pi)  
        omega2 = omega2 - pi  
   
        TEST_OMEGA12: IF (omega1<=omega2) THEN  
           ! on est dans la meme journee locale  
           IF (omega2<=-omega .OR. omega1>=omega .OR. omega<1E-5) THEN  
              ! nuit  
              IF (present(frac)) frac(i) = 0.0  
              pmu0(i) = 0.0  
           ELSE  
              ! jour + nuit / jour  
              omegadeb = max(-omega, omega1)  
              omegafin = min(omega, omega2)  
              IF (present(frac)) frac(i) = (omegafin-omegadeb)/(omega2-omega1)  
              pmu0(i) = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &  
                   omegafin)-sin(omegadeb))/(omegafin-omegadeb)  
           END IF  
        ELSE TEST_OMEGA12  
           !---omega1 GT omega2 -- a cheval sur deux journees  
           !-------------------entre omega1 et pi  
           IF (omega1>=omega) THEN !--nuit  
              zfrac1 = 0.0  
              z1_mu = 0.0  
           ELSE !--jour+nuit  
              omegadeb = max(-omega, omega1)  
              omegafin = omega  
              zfrac1 = omegafin - omegadeb  
              z1_mu = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &  
                   omegafin)-sin(omegadeb))/(omegafin-omegadeb)  
           END IF  
           !---------------------entre -pi et omega2  
           IF (omega2<=-omega) THEN !--nuit  
              zfrac2 = 0.0  
              z2_mu = 0.0  
           ELSE !--jour+nuit  
              omegadeb = -omega  
              omegafin = min(omega, omega2)  
              zfrac2 = omegafin - omegadeb  
              z2_mu = sin(latr)*sin(lat_sun) + cos(latr)*cos(lat_sun)*(sin( &  
                   omegafin)-sin(omegadeb))/(omegafin-omegadeb)  
   
           END IF  
           !-----------------------moyenne  
           IF (present(frac)) frac(i) = (zfrac1+zfrac2)/ &  
                (omega2+deux_pi-omega1)  
           pmu0(i) = (zfrac1*z1_mu+zfrac2*z2_mu)/max(zfrac1+zfrac2, 1.E-10)  
        END IF TEST_OMEGA12  
     END DO  
   
   END SUBROUTINE zenang  
   
61  END MODULE orbite_m  END MODULE orbite_m
Legend:
Removed from v.36  
changed lines
  Added in v.118
  ViewVC Help
Powered by ViewVC 1.1.21