trunk/phylmd/albedo.f


! $Header: /home/cvsroot/LMDZ4/libf/phylmd/albedo.F,v 1.2 2005/02/07 15:00:52
! fairhead Exp $


SUBROUTINE alboc(rjour, rlat, albedo)
  USE dimens_m
  USE dimphy
  USE yomcst
  USE orbite_m, ONLY: orbite
  IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
  ! Date: le 16 mars 1995
  ! Objet: Calculer l'albedo sur l'ocean
  ! Methode: Integrer numeriquement l'albedo pendant une journee

  ! Arguments;
  ! rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)
  ! rlat (in,R)   : latitude en degre
  ! albedo (out,R): albedo obtenu (de 0 a 1)
  ! ======================================================================

  REAL fmagic ! un facteur magique pour regler l'albedo
  ! cc      PARAMETER (fmagic=0.7)
  ! ccIM => a remplacer
  ! PARAMETER (fmagic=1.32)
  PARAMETER (fmagic=1.0)
  ! PARAMETER (fmagic=0.7)
  INTEGER npts ! il controle la precision de l'integration
  PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes

  REAL rlat(klon), rjour, albedo(klon)
  REAL zdist, zlonsun, zpi, zdeclin
  REAL rmu, alb, srmu, salb, fauxo, aa, bb
  INTEGER i, k
  ! ccIM
  LOGICAL ancien_albedo
  PARAMETER (ancien_albedo=.FALSE.)
  ! SAVE albedo

  IF (ancien_albedo) THEN

    zpi = 4.*atan(1.)

    ! Calculer la longitude vraie de l'orbite terrestre:
    CALL orbite(rjour, zlonsun, zdist)

    ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
    zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))

    DO i = 1, klon
      aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
      bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)

      ! 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

      ! Faire l'integration numerique de midi a minuit (le facteur 2
      ! 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
      END DO
      IF (srmu/=0.0) THEN
        albedo(i) = salb/srmu*fmagic
      ELSE ! nuit polaire (on peut prendre une valeur quelconque)
        albedo(i) = fmagic
      END IF
    END DO

    ! nouvel albedo

  ELSE

    zpi = 4.*atan(1.)

    ! Calculer la longitude vraie de l'orbite terrestre:
    CALL orbite(rjour, zlonsun, zdist)

    ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
    zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))

    DO i = 1, klon
      aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
      bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)

      ! Midi local (angle du temps = 0.0):
      rmu = aa + bb*cos(0.0)
      rmu = max(0.0, rmu)
      ! IM cf. PB  alb = 0.058/(rmu + 0.30)
      ! alb = 0.058/(rmu + 0.30) * 1.5
      alb = 0.058/(rmu+0.30)*1.2
      ! alb = 0.058/(rmu + 0.30) * 1.3
      srmu = rmu
      salb = alb*rmu

      ! Faire l'integration numerique de midi a minuit (le facteur 2
      ! 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)
        ! IM cf. PB      alb = 0.058/(rmu + 0.30)
        ! alb = 0.058/(rmu + 0.30) * 1.5
        alb = 0.058/(rmu+0.30)*1.2
        ! alb = 0.058/(rmu + 0.30) * 1.3
        srmu = srmu + rmu*2.0
        salb = salb + alb*rmu*2.0
      END DO
      IF (srmu/=0.0) THEN
        albedo(i) = salb/srmu*fmagic
      ELSE ! nuit polaire (on peut prendre une valeur quelconque)
        albedo(i) = fmagic
      END IF
    END DO
  END IF
  RETURN
END SUBROUTINE alboc
! =====================================================================
SUBROUTINE alboc_cd(rmu0, albedo)
  USE dimens_m
  USE dimphy
  IMPLICIT NONE
  ! ======================================================================
  ! Auteur(s): Z.X. Li (LMD/CNRS)
  ! date: 19940624
  ! Calculer l'albedo sur l'ocean en fonction de l'angle zenithal moyen
  ! Formule due a Larson and Barkstrom (1977) Proc. of the symposium
  ! on radiation in the atmosphere, 19-28 August 1976, science Press,
  ! 1977 pp 451-453, ou These de 3eme cycle de Sylvie Joussaume.

  ! Arguments
  ! rmu0    (in): cosinus de l'angle solaire zenithal
  ! albedo (out): albedo de surface de l'ocean
  ! ======================================================================
  REAL rmu0(klon), albedo(klon)

  REAL fmagic ! un facteur magique pour regler l'albedo
  ! cc      PARAMETER (fmagic=0.7)
  ! ccIM => a remplacer
  ! PARAMETER (fmagic=1.32)
  PARAMETER (fmagic=1.0)
  ! PARAMETER (fmagic=0.7)

  REAL fauxo
  INTEGER i
  ! ccIM
  LOGICAL ancien_albedo
  PARAMETER (ancien_albedo=.FALSE.)
  ! SAVE albedo

  IF (ancien_albedo) THEN

    DO i = 1, klon

      rmu0(i) = max(rmu0(i), 0.0)

      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)
    END DO

    ! nouvel albedo

  ELSE

    DO i = 1, klon
      rmu0(i) = max(rmu0(i), 0.0)
      ! IM: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)
    END DO

  END IF

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