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

Diff of /trunk/phylmd/alboc.f

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

trunk/libf/phylmd/albedo.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/phylmd/albedo.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/albedo.F,v 1.2 2005/02/07 15:00:52 fairhead Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/albedo.F,v 1.2 2005/02/07 15:00:52
3  !  ! fairhead Exp $
4  c  
5  c  
6        SUBROUTINE alboc(rjour,rlat,albedo)  
7        use dimens_m  SUBROUTINE alboc(rjour, rlat, albedo)
8        use dimphy    USE dimens_m
9        use YOMCST    USE dimphy
10        use orbite_m, only: orbite    USE yomcst
11        IMPLICIT none    USE orbite_m, ONLY: orbite
12  c======================================================================    IMPLICIT NONE
13  c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)    ! ======================================================================
14  c Date: le 16 mars 1995    ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
15  c Objet: Calculer l'albedo sur l'ocean    ! Date: le 16 mars 1995
16  c Methode: Integrer numeriquement l'albedo pendant une journee    ! Objet: Calculer l'albedo sur l'ocean
17  c    ! Methode: Integrer numeriquement l'albedo pendant une journee
18  c Arguments;  
19  c rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)    ! Arguments;
20  c rlat (in,R)   : latitude en degre    ! rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)
21  c albedo (out,R): albedo obtenu (de 0 a 1)    ! rlat (in,R)   : latitude en degre
22  c======================================================================    ! albedo (out,R): albedo obtenu (de 0 a 1)
23  c    ! ======================================================================
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    REAL rlat(klon), rjour, albedo(klon)
35        REAL rmu,alb, srmu, salb, fauxo, aa, bb    REAL zdist, zlonsun, zpi, zdeclin
36        INTEGER i, k    REAL rmu, alb, srmu, salb, fauxo, aa, bb
37  cccIM    INTEGER i, k
38        LOGICAL ancien_albedo    ! ccIM
39        PARAMETER(ancien_albedo=.FALSE.)    LOGICAL ancien_albedo
40  c     SAVE albedo    PARAMETER (ancien_albedo=.FALSE.)
41  c    ! SAVE albedo
42        IF ( ancien_albedo ) THEN  
43  c    IF (ancien_albedo) THEN
44        zpi = 4. * ATAN(1.)  
45  c      zpi = 4.*atan(1.)
46  c Calculer la longitude vraie de l'orbite terrestre:  
47        CALL orbite(rjour,zlonsun,zdist)      ! Calculer la longitude vraie de l'orbite terrestre:
48  c      CALL orbite(rjour, zlonsun, zdist)
49  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):  
50        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))      ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
51  c      zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
52        DO 999 i=1,klon  
53        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)      DO i = 1, klon
54        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)        aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
55  c        bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
56  c Midi local (angle du temps = 0.0):  
57        rmu = aa + bb * COS(0.0)        ! Midi local (angle du temps = 0.0):
58        rmu = MAX(0.0, rmu)        rmu = aa + bb*cos(0.0)
59        fauxo = (1.47-ACOS(rmu))/.15        rmu = max(0.0, rmu)
60        alb = 0.03+0.630/(1.+fauxo*fauxo)        fauxo = (1.47-acos(rmu))/.15
61          alb = 0.03 + 0.630/(1.+fauxo*fauxo)
62        srmu = rmu        srmu = rmu
63        salb = alb * rmu        salb = alb*rmu
64  c  
65  c Faire l'integration numerique de midi a minuit (le facteur 2        ! Faire l'integration numerique de midi a minuit (le facteur 2
66  c prend en compte l'autre moitie de la journee):        ! prend en compte l'autre moitie de la journee):
67        DO k = 1, npts        DO k = 1, npts
68           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)          rmu = aa + bb*cos(float(k)/float(npts)*zpi)
69           rmu = MAX(0.0, rmu)          rmu = max(0.0, rmu)
70           fauxo = (1.47-ACOS(rmu))/.15          fauxo = (1.47-acos(rmu))/.15
71           alb = 0.03+0.630/(1.+fauxo*fauxo)          alb = 0.03 + 0.630/(1.+fauxo*fauxo)
72           srmu = srmu + rmu * 2.0          srmu = srmu + rmu*2.0
73           salb = salb + alb*rmu * 2.0          salb = salb + alb*rmu*2.0
74        ENDDO        END DO
75        IF (srmu .NE. 0.0) THEN        IF (srmu/=0.0) THEN
76           albedo(i) = salb / srmu * fmagic          albedo(i) = salb/srmu*fmagic
77        ELSE ! nuit polaire (on peut prendre une valeur quelconque)        ELSE ! nuit polaire (on peut prendre une valeur quelconque)
78           albedo(i) = fmagic          albedo(i) = fmagic
79        ENDIF        END IF
80    999 CONTINUE      END DO
81  c  
82  c nouvel albedo      ! nouvel albedo
83  c  
84        ELSE    ELSE
85  c  
86        zpi = 4. * ATAN(1.)      zpi = 4.*atan(1.)
87  c  
88  c Calculer la longitude vraie de l'orbite terrestre:      ! Calculer la longitude vraie de l'orbite terrestre:
89        CALL orbite(rjour,zlonsun,zdist)      CALL orbite(rjour, zlonsun, zdist)
90  c  
91  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):      ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
92        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))      zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
93  c  
94        DO 1999 i=1,klon      DO i = 1, klon
95        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)        aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
96        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)        bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
97  c  
98  c Midi local (angle du temps = 0.0):        ! Midi local (angle du temps = 0.0):
99        rmu = aa + bb * COS(0.0)        rmu = aa + bb*cos(0.0)
100        rmu = MAX(0.0, rmu)        rmu = max(0.0, rmu)
101  cIM cf. PB  alb = 0.058/(rmu + 0.30)        ! IM cf. PB  alb = 0.058/(rmu + 0.30)
102  c     alb = 0.058/(rmu + 0.30) * 1.5        ! alb = 0.058/(rmu + 0.30) * 1.5
103        alb = 0.058/(rmu + 0.30) * 1.2        alb = 0.058/(rmu+0.30)*1.2
104  c     alb = 0.058/(rmu + 0.30) * 1.3        ! alb = 0.058/(rmu + 0.30) * 1.3
105        srmu = rmu        srmu = rmu
106        salb = alb * rmu        salb = alb*rmu
107  c  
108  c Faire l'integration numerique de midi a minuit (le facteur 2        ! Faire l'integration numerique de midi a minuit (le facteur 2
109  c prend en compte l'autre moitie de la journee):        ! prend en compte l'autre moitie de la journee):
110        DO k = 1, npts        DO k = 1, npts
111           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)          rmu = aa + bb*cos(float(k)/float(npts)*zpi)
112           rmu = MAX(0.0, rmu)          rmu = max(0.0, rmu)
113  cIM cf. PB      alb = 0.058/(rmu + 0.30)          ! IM cf. PB      alb = 0.058/(rmu + 0.30)
114  c        alb = 0.058/(rmu + 0.30) * 1.5          ! alb = 0.058/(rmu + 0.30) * 1.5
115           alb = 0.058/(rmu + 0.30) * 1.2          alb = 0.058/(rmu+0.30)*1.2
116  c        alb = 0.058/(rmu + 0.30) * 1.3          ! alb = 0.058/(rmu + 0.30) * 1.3
117           srmu = srmu + rmu * 2.0          srmu = srmu + rmu*2.0
118           salb = salb + alb*rmu * 2.0          salb = salb + alb*rmu*2.0
119        ENDDO        END DO
120        IF (srmu .NE. 0.0) THEN        IF (srmu/=0.0) THEN
121           albedo(i) = salb / srmu * fmagic          albedo(i) = salb/srmu*fmagic
122        ELSE ! nuit polaire (on peut prendre une valeur quelconque)        ELSE ! nuit polaire (on peut prendre une valeur quelconque)
123           albedo(i) = fmagic          albedo(i) = fmagic
124        ENDIF        END IF
125  1999  CONTINUE      END DO
126        ENDIF    END IF
127        RETURN    RETURN
128        END  END SUBROUTINE alboc
129  c=====================================================================  ! =====================================================================
130        SUBROUTINE alboc_cd(rmu0,albedo)  SUBROUTINE alboc_cd(rmu0, albedo)
131        use dimens_m    USE dimens_m
132        use dimphy    USE dimphy
133        IMPLICIT none    IMPLICIT NONE
134  c======================================================================    ! ======================================================================
135  c Auteur(s): Z.X. Li (LMD/CNRS)    ! Auteur(s): Z.X. Li (LMD/CNRS)
136  c date: 19940624    ! date: 19940624
137  c Calculer l'albedo sur l'ocean en fonction de l'angle zenithal moyen    ! Calculer l'albedo sur l'ocean en fonction de l'angle zenithal moyen
138  c Formule due a Larson and Barkstrom (1977) Proc. of the symposium    ! Formule due a Larson and Barkstrom (1977) Proc. of the symposium
139  C on radiation in the atmosphere, 19-28 August 1976, science Press,    ! on radiation in the atmosphere, 19-28 August 1976, science Press,
140  C 1977 pp 451-453, ou These de 3eme cycle de Sylvie Joussaume.    ! 1977 pp 451-453, ou These de 3eme cycle de Sylvie Joussaume.
141  c  
142  c Arguments    ! Arguments
143  c rmu0    (in): cosinus de l'angle solaire zenithal    ! rmu0    (in): cosinus de l'angle solaire zenithal
144  c albedo (out): albedo de surface de l'ocean    ! albedo (out): albedo de surface de l'ocean
145  c======================================================================    ! ======================================================================
146        REAL rmu0(klon), albedo(klon)    REAL rmu0(klon), albedo(klon)
147  c  
148        REAL fmagic ! un facteur magique pour regler l'albedo    REAL fmagic ! un facteur magique pour regler l'albedo
149  ccc      PARAMETER (fmagic=0.7)    ! cc      PARAMETER (fmagic=0.7)
150  cccIM => a remplacer      ! ccIM => a remplacer
151  c       PARAMETER (fmagic=1.32)    ! PARAMETER (fmagic=1.32)
152          PARAMETER (fmagic=1.0)    PARAMETER (fmagic=1.0)
153  c       PARAMETER (fmagic=0.7)    ! PARAMETER (fmagic=0.7)
154  c  
155        REAL fauxo    REAL fauxo
156        INTEGER i    INTEGER i
157  cccIM    ! ccIM
158        LOGICAL ancien_albedo    LOGICAL ancien_albedo
159        PARAMETER(ancien_albedo=.FALSE.)    PARAMETER (ancien_albedo=.FALSE.)
160  c     SAVE albedo    ! SAVE albedo
161  c  
162        IF ( ancien_albedo ) THEN    IF (ancien_albedo) THEN
163  c  
164        DO i = 1, klon      DO i = 1, klon
165  c  
166           rmu0(i) = MAX(rmu0(i),0.0)        rmu0(i) = max(rmu0(i), 0.0)
167  c  
168           fauxo = ( 1.47 - ACOS( rmu0(i) ) )/0.15        fauxo = (1.47-acos(rmu0(i)))/0.15
169           albedo(i) = fmagic*( .03 + .630/( 1. + fauxo*fauxo))        albedo(i) = fmagic*(.03+.630/(1.+fauxo*fauxo))
170           albedo(i) = MAX(MIN(albedo(i),0.60),0.04)        albedo(i) = max(min(albedo(i),0.60), 0.04)
171        ENDDO      END DO
172  c  
173  c nouvel albedo      ! nouvel albedo
174  c  
175        ELSE    ELSE
176  c  
177        DO i = 1, klon      DO i = 1, klon
178           rmu0(i) = MAX(rmu0(i),0.0)        rmu0(i) = max(rmu0(i), 0.0)
179  cIM:orig albedo(i) = 0.058/(rmu0(i) + 0.30)        ! IM:orig albedo(i) = 0.058/(rmu0(i) + 0.30)
180           albedo(i) = fmagic * 0.058/(rmu0(i) + 0.30)        albedo(i) = fmagic*0.058/(rmu0(i)+0.30)
181           albedo(i) = MAX(MIN(albedo(i),0.60),0.04)        albedo(i) = max(min(albedo(i),0.60), 0.04)
182        ENDDO      END DO
183  c  
184        ENDIF    END IF
185  c  
186        RETURN    RETURN
187        END  END SUBROUTINE alboc_cd
188  c========================================================================  ! ========================================================================
Legend:
Removed from v.3  
changed lines
  Added in v.82
  ViewVC Help
Powered by ViewVC 1.1.21