source: LMDZ4/branches/pre_V3/libf/phylmd/albedo.F @ 1365

Last change on this file since 1365 was 584, checked in by Laurent Fairhead, 20 years ago

Pour se mettre en conformite avec les simulation IPCC
LF

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