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