libf/phylmd/albedo.f

!
! $Header: /home/cvsroot/LMDZ4/libf/phylmd/albedo.F,v 1.2 2005/02/07 15:00:52 fairhead Exp $
!
c
c
      SUBROUTINE alboc(rjour,rlat,albedo)
      use dimens_m
      use dimphy
      use yomcst
      use orbite_m, only: orbite
      IMPLICIT none
c======================================================================
c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
c Date: le 16 mars 1995
c Objet: Calculer l'albedo sur l'ocean
c Methode: Integrer numeriquement l'albedo pendant une journee
c
c Arguments;
c rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)
c rlat (in,R)   : latitude en degre
c albedo (out,R): albedo obtenu (de 0 a 1)
c======================================================================
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)
      INTEGER npts ! il controle la precision de l'integration
      PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes
c
      REAL rlat(klon), rjour, albedo(klon)
      REAL zdist, zlonsun, zpi, zdeclin
      REAL rmu,alb, srmu, salb, fauxo, aa, bb
      INTEGER i, k
cccIM
      LOGICAL ancien_albedo
      PARAMETER(ancien_albedo=.FALSE.) 
c     SAVE albedo
c
      IF ( ancien_albedo ) THEN
c
      zpi = 4. * ATAN(1.)
c
c Calculer la longitude vraie de l'orbite terrestre:
      CALL orbite(rjour,zlonsun,zdist)
c
c Calculer la declinaison du soleil (qui varie entre + et - R_incl):
      zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))
c
      DO 999 i=1,klon
      aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)
      bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)
c
c Midi local (angle du temps = 0.0):
      rmu = aa + bb * COS(0.0)
      rmu = MAX(0.0, rmu)
      fauxo = (1.47-ACOS(rmu))/.15
      alb = 0.03+0.630/(1.+fauxo*fauxo)
      srmu = rmu
      salb = alb * rmu
c
c Faire l'integration numerique de midi a minuit (le facteur 2
c prend en compte l'autre moitie de la journee):
      DO k = 1, npts
         rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)
         rmu = MAX(0.0, rmu)
         fauxo = (1.47-ACOS(rmu))/.15
         alb = 0.03+0.630/(1.+fauxo*fauxo)
         srmu = srmu + rmu * 2.0
         salb = salb + alb*rmu * 2.0
      ENDDO
      IF (srmu .NE. 0.0) THEN
         albedo(i) = salb / srmu * fmagic
      ELSE ! nuit polaire (on peut prendre une valeur quelconque)
         albedo(i) = fmagic
      ENDIF
  999 CONTINUE
c
c nouvel albedo 
c
      ELSE
c
      zpi = 4. * ATAN(1.)
c
c Calculer la longitude vraie de l'orbite terrestre:
      CALL orbite(rjour,zlonsun,zdist)
c
c Calculer la declinaison du soleil (qui varie entre + et - R_incl):
      zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))
c
      DO 1999 i=1,klon
      aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)
      bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)
c
c Midi local (angle du temps = 0.0):
      rmu = aa + bb * COS(0.0)
      rmu = MAX(0.0, rmu)
cIM cf. PB  alb = 0.058/(rmu + 0.30)
c     alb = 0.058/(rmu + 0.30) * 1.5
      alb = 0.058/(rmu + 0.30) * 1.2
c     alb = 0.058/(rmu + 0.30) * 1.3
      srmu = rmu
      salb = alb * rmu
c
c Faire l'integration numerique de midi a minuit (le facteur 2
c prend en compte l'autre moitie de la journee):
      DO k = 1, npts
         rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)
         rmu = MAX(0.0, rmu)
cIM cf. PB      alb = 0.058/(rmu + 0.30)
c        alb = 0.058/(rmu + 0.30) * 1.5
         alb = 0.058/(rmu + 0.30) * 1.2
c        alb = 0.058/(rmu + 0.30) * 1.3
         srmu = srmu + rmu * 2.0
         salb = salb + alb*rmu * 2.0
      ENDDO
      IF (srmu .NE. 0.0) THEN
         albedo(i) = salb / srmu * fmagic
      ELSE ! nuit polaire (on peut prendre une valeur quelconque)
         albedo(i) = fmagic
      ENDIF
1999  CONTINUE
      ENDIF
      RETURN
      END
c=====================================================================
      SUBROUTINE alboc_cd(rmu0,albedo)
      use dimens_m
      use dimphy
      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========================================================================
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(zlonsunzpi/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 + albrmu 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(zlonsunzpi/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 + albrmu 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. + fauxofauxo))
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========================================================================