/[lmdze]/trunk/phylmd/Interface_surf/alboc.f
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/alboc.f

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

trunk/libf/phylmd/albedo.f revision 38 by guez, Thu Jan 6 17:52:19 2011 UTC trunk/phylmd/Interface_surf/alboc.f revision 130 by guez, Tue Feb 24 15:43:51 2015 UTC
# Line 1  Line 1 
1  !  module alboc_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/albedo.F,v 1.2 2005/02/07 15:00:52 fairhead Exp $  
3  !    IMPLICIT NONE
4  c  
5  c  contains
6        SUBROUTINE alboc(rjour,rlat,albedo)  
7        use dimens_m    SUBROUTINE alboc(jour, rlat, albedo)
8        use dimphy  
9        use yomcst      ! From LMDZ4/libf/phylmd/albedo.F, version 1.2 2005/02/07 15:00:52
10        use orbite_m, only: orbite  
11        IMPLICIT none      ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
12  c======================================================================      ! Date: le 16 mars 1995
13  c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)      ! Objet: Calculer l'albedo sur l'ocean
14  c Date: le 16 mars 1995      ! Methode: Integrer numeriquement l'albedo pendant une journee
15  c Objet: Calculer l'albedo sur l'ocean  
16  c Methode: Integrer numeriquement l'albedo pendant une journee      USE dimphy, only: klon
17  c      USE yomcst, only: r_incl
18  c Arguments;      USE orbite_m, ONLY: orbite
19  c rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)  
20  c rlat (in,R)   : latitude en degre      ! Arguments;
21  c albedo (out,R): albedo obtenu (de 0 a 1)      ! jour (in) : jour dans l'annee (a compter du 1 janvier)
22  c======================================================================      ! rlat (in, R) : latitude en degre
23  c      ! albedo (out, R): albedo obtenu (de 0 a 1)
24        REAL fmagic ! un facteur magique pour regler l'albedo  
25  ccc      PARAMETER (fmagic=0.7)      REAL fmagic ! un facteur magique pour regler l'albedo
26  cccIM => a remplacer        ! cc PARAMETER (fmagic=0.7)
27  c       PARAMETER (fmagic=1.32)      ! ccIM => a remplacer
28          PARAMETER (fmagic=1.0)      ! PARAMETER (fmagic=1.32)
29  c       PARAMETER (fmagic=0.7)      PARAMETER (fmagic=1.0)
30        INTEGER npts ! il controle la precision de l'integration      ! PARAMETER (fmagic=0.7)
31        PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes      INTEGER npts ! il controle la precision de l'integration
32  c      PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes
33        REAL rlat(klon), rjour, albedo(klon)  
34        REAL zdist, zlonsun, zpi, zdeclin      integer jour
35        REAL rmu,alb, srmu, salb, fauxo, aa, bb      REAL rlat(klon), albedo(klon)
36        INTEGER i, k      REAL zdist, zlonsun, zpi, zdeclin
37  cccIM      REAL rmu, alb, srmu, salb, fauxo, aa, bb
38        LOGICAL ancien_albedo      INTEGER i, k
39        PARAMETER(ancien_albedo=.FALSE.)      ! ccIM
40  c     SAVE albedo      LOGICAL ancien_albedo
41  c      PARAMETER (ancien_albedo=.FALSE.)
42        IF ( ancien_albedo ) THEN      ! SAVE albedo
43  c  
44        zpi = 4. * ATAN(1.)      IF (ancien_albedo) THEN
45  c  
46  c Calculer la longitude vraie de l'orbite terrestre:         zpi = 4.*atan(1.)
47        CALL orbite(rjour,zlonsun,zdist)  
48  c         ! Calculer la longitude vraie de l'orbite terrestre:
49  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):         CALL orbite(real(jour), zlonsun, zdist)
50        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))  
51  c         ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
52        DO 999 i=1,klon         zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
53        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)  
54        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)         DO i = 1, klon
55  c            aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
56  c Midi local (angle du temps = 0.0):            bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
57        rmu = aa + bb * COS(0.0)  
58        rmu = MAX(0.0, rmu)            ! Midi local (angle du temps = 0.0):
59        fauxo = (1.47-ACOS(rmu))/.15            rmu = aa + bb*cos(0.0)
60        alb = 0.03+0.630/(1.+fauxo*fauxo)            rmu = max(0.0, rmu)
61        srmu = rmu            fauxo = (1.47-acos(rmu))/.15
62        salb = alb * rmu            alb = 0.03 + 0.630/(1.+fauxo*fauxo)
63  c            srmu = rmu
64  c Faire l'integration numerique de midi a minuit (le facteur 2            salb = alb*rmu
65  c prend en compte l'autre moitie de la journee):  
66        DO k = 1, npts            ! Faire l'integration numerique de midi a minuit (le facteur 2
67           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)            ! prend en compte l'autre moitie de la journee):
68           rmu = MAX(0.0, rmu)            DO k = 1, npts
69           fauxo = (1.47-ACOS(rmu))/.15               rmu = aa + bb*cos(float(k)/float(npts)*zpi)
70           alb = 0.03+0.630/(1.+fauxo*fauxo)               rmu = max(0.0, rmu)
71           srmu = srmu + rmu * 2.0               fauxo = (1.47-acos(rmu))/.15
72           salb = salb + alb*rmu * 2.0               alb = 0.03 + 0.630/(1.+fauxo*fauxo)
73        ENDDO               srmu = srmu + rmu*2.0
74        IF (srmu .NE. 0.0) THEN               salb = salb + alb*rmu*2.0
75           albedo(i) = salb / srmu * fmagic            END DO
76        ELSE ! nuit polaire (on peut prendre une valeur quelconque)            IF (srmu/=0.0) THEN
77           albedo(i) = fmagic               albedo(i) = salb/srmu*fmagic
78        ENDIF            ELSE ! nuit polaire (on peut prendre une valeur quelconque)
79    999 CONTINUE               albedo(i) = fmagic
80  c            END IF
81  c nouvel albedo         END DO
82  c  
83        ELSE         ! nouvel albedo
84  c  
85        zpi = 4. * ATAN(1.)      ELSE
86  c  
87  c Calculer la longitude vraie de l'orbite terrestre:         zpi = 4.*atan(1.)
88        CALL orbite(rjour,zlonsun,zdist)  
89  c         ! Calculer la longitude vraie de l'orbite terrestre:
90  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):         CALL orbite(real(jour), zlonsun, zdist)
91        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))  
92  c         ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
93        DO 1999 i=1,klon         zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
94        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)  
95        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)         DO i = 1, klon
96  c            aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
97  c Midi local (angle du temps = 0.0):            bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
98        rmu = aa + bb * COS(0.0)  
99        rmu = MAX(0.0, rmu)            ! Midi local (angle du temps = 0.0):
100  cIM cf. PB  alb = 0.058/(rmu + 0.30)            rmu = aa + bb*cos(0.0)
101  c     alb = 0.058/(rmu + 0.30) * 1.5            rmu = max(0.0, rmu)
102        alb = 0.058/(rmu + 0.30) * 1.2            ! IM cf. PB alb = 0.058/(rmu + 0.30)
103  c     alb = 0.058/(rmu + 0.30) * 1.3            ! alb = 0.058/(rmu + 0.30) * 1.5
104        srmu = rmu            alb = 0.058/(rmu+0.30)*1.2
105        salb = alb * rmu            ! alb = 0.058/(rmu + 0.30) * 1.3
106  c            srmu = rmu
107  c Faire l'integration numerique de midi a minuit (le facteur 2            salb = alb*rmu
108  c prend en compte l'autre moitie de la journee):  
109        DO k = 1, npts            ! Faire l'integration numerique de midi a minuit (le facteur 2
110           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)            ! prend en compte l'autre moitie de la journee):
111           rmu = MAX(0.0, rmu)            DO k = 1, npts
112  cIM cf. PB      alb = 0.058/(rmu + 0.30)               rmu = aa + bb*cos(float(k)/float(npts)*zpi)
113  c        alb = 0.058/(rmu + 0.30) * 1.5               rmu = max(0.0, rmu)
114           alb = 0.058/(rmu + 0.30) * 1.2               ! IM cf. PB alb = 0.058/(rmu + 0.30)
115  c        alb = 0.058/(rmu + 0.30) * 1.3               ! alb = 0.058/(rmu + 0.30) * 1.5
116           srmu = srmu + rmu * 2.0               alb = 0.058/(rmu+0.30)*1.2
117           salb = salb + alb*rmu * 2.0               ! alb = 0.058/(rmu + 0.30) * 1.3
118        ENDDO               srmu = srmu + rmu*2.0
119        IF (srmu .NE. 0.0) THEN               salb = salb + alb*rmu*2.0
120           albedo(i) = salb / srmu * fmagic            END DO
121        ELSE ! nuit polaire (on peut prendre une valeur quelconque)            IF (srmu/=0.0) THEN
122           albedo(i) = fmagic               albedo(i) = salb/srmu*fmagic
123        ENDIF            ELSE ! nuit polaire (on peut prendre une valeur quelconque)
124  1999  CONTINUE               albedo(i) = fmagic
125        ENDIF            END IF
126        RETURN         END DO
127        END      END IF
128  c=====================================================================  
129        SUBROUTINE alboc_cd(rmu0,albedo)    END SUBROUTINE alboc
130        use dimens_m  
131        use dimphy  end module alboc_m
       IMPLICIT none  
 c======================================================================  
 c Auteur(s): Z.X. Li (LMD/CNRS)  
 c date: 19940624  
 c Calculer l'albedo sur l'ocean en fonction de l'angle zenithal moyen  
 c Formule due a Larson and Barkstrom (1977) Proc. of the symposium  
 C on radiation in the atmosphere, 19-28 August 1976, science Press,  
 C 1977 pp 451-453, ou These de 3eme cycle de Sylvie Joussaume.  
 c  
 c Arguments  
 c rmu0    (in): cosinus de l'angle solaire zenithal  
 c albedo (out): albedo de surface de l'ocean  
 c======================================================================  
       REAL rmu0(klon), albedo(klon)  
 c  
       REAL fmagic ! un facteur magique pour regler l'albedo  
 ccc      PARAMETER (fmagic=0.7)  
 cccIM => a remplacer    
 c       PARAMETER (fmagic=1.32)  
         PARAMETER (fmagic=1.0)  
 c       PARAMETER (fmagic=0.7)  
 c  
       REAL fauxo  
       INTEGER i  
 cccIM  
       LOGICAL ancien_albedo  
       PARAMETER(ancien_albedo=.FALSE.)  
 c     SAVE albedo  
 c  
       IF ( ancien_albedo ) THEN  
 c  
       DO i = 1, klon  
 c  
          rmu0(i) = MAX(rmu0(i),0.0)  
 c  
          fauxo = ( 1.47 - ACOS( rmu0(i) ) )/0.15  
          albedo(i) = fmagic*( .03 + .630/( 1. + fauxo*fauxo))  
          albedo(i) = MAX(MIN(albedo(i),0.60),0.04)  
       ENDDO  
 c  
 c nouvel albedo  
 c  
       ELSE  
 c  
       DO i = 1, klon  
          rmu0(i) = MAX(rmu0(i),0.0)  
 cIM:orig albedo(i) = 0.058/(rmu0(i) + 0.30)  
          albedo(i) = fmagic * 0.058/(rmu0(i) + 0.30)  
          albedo(i) = MAX(MIN(albedo(i),0.60),0.04)  
       ENDDO  
 c  
       ENDIF  
 c  
       RETURN  
       END  
 c========================================================================  
Legend:
Removed from v.38  
changed lines
  Added in v.130
  ViewVC Help
Powered by ViewVC 1.1.21