phylmd/Interface_surf/alboc.f

module alboc_m

  IMPLICIT NONE

contains

  SUBROUTINE alboc(jour, rlat, albedo)

    ! From LMDZ4/libf/phylmd/albedo.F, version 1.2 2005/02/07 15:00:52

    ! 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

    USE dimphy, only: klon
    USE yomcst, only: r_incl
    USE orbite_m, ONLY: orbite

    ! Arguments;
    ! jour (in) : 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

    integer jour
    REAL rlat(klon), 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(real(jour), 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(real(jour), 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

  END SUBROUTINE alboc

end module alboc_m
1	guez	125	module alboc_m
2
3	guez	81	IMPLICIT NONE
4
5	guez	125	contains
6	guez	81
7	guez	130	SUBROUTINE alboc(jour, rlat, albedo)
8	guez	81
9	guez	125	! From LMDZ4/libf/phylmd/albedo.F, version 1.2 2005/02/07 15:00:52
10	guez	81
11	guez	125	! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
12			! Date: le 16 mars 1995
13			! Objet: Calculer l'albedo sur l'ocean
14			! Methode: Integrer numeriquement l'albedo pendant une journee
15	guez	81
16	guez	125	USE dimphy, only: klon
17			USE yomcst, only: r_incl
18			USE orbite_m, ONLY: orbite
19	guez	81
20	guez	125	! Arguments;
21	guez	130	! jour (in) : jour dans l'annee (a compter du 1 janvier)
22	guez	125	! rlat (in, R) : latitude en degre
23			! albedo (out, R): albedo obtenu (de 0 a 1)
24	guez	81
25	guez	125	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	guez	81
34	guez	130	integer jour
35			REAL rlat(klon), albedo(klon)
36	guez	125	REAL zdist, zlonsun, zpi, zdeclin
37			REAL rmu, alb, srmu, salb, fauxo, aa, bb
38			INTEGER i, k
39			! ccIM
40			LOGICAL ancien_albedo
41			PARAMETER (ancien_albedo=.FALSE.)
42			! SAVE albedo
43	guez	81
44	guez	125	IF (ancien_albedo) THEN
45	guez	81
46	guez	125	zpi = 4.*atan(1.)
47	guez	81
48	guez	125	! Calculer la longitude vraie de l'orbite terrestre:
49	guez	130	CALL orbite(real(jour), zlonsun, zdist)
50	guez	81
51	guez	125	! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
52			zdeclin = asin(sin(zlonsunzpi/180.0)sin(r_incl*zpi/180.0))
53	guez	81
54	guez	125	DO i = 1, klon
55			aa = sin(rlat(i)zpi/180.0)sin(zdeclin)
56			bb = cos(rlat(i)zpi/180.0)cos(zdeclin)
57	guez	81
58	guez	125	! Midi local (angle du temps = 0.0):
59			rmu = aa + bb*cos(0.0)
60			rmu = max(0.0, rmu)
61			fauxo = (1.47-acos(rmu))/.15
62			alb = 0.03 + 0.630/(1.+fauxo*fauxo)
63			srmu = rmu
64			salb = alb*rmu
65	guez	81
66	guez	125	! Faire l'integration numerique de midi a minuit (le facteur 2
67			! prend en compte l'autre moitie de la journee):
68			DO k = 1, npts
69			rmu = aa + bbcos(float(k)/float(npts)zpi)
70			rmu = max(0.0, rmu)
71			fauxo = (1.47-acos(rmu))/.15
72			alb = 0.03 + 0.630/(1.+fauxo*fauxo)
73			srmu = srmu + rmu*2.0
74			salb = salb + albrmu2.0
75			END DO
76			IF (srmu/=0.0) THEN
77			albedo(i) = salb/srmu*fmagic
78			ELSE ! nuit polaire (on peut prendre une valeur quelconque)
79			albedo(i) = fmagic
80			END IF
81			END DO
82	guez	81
83	guez	125	! nouvel albedo
84	guez	81
85	guez	125	ELSE
86	guez	81
87	guez	125	zpi = 4.*atan(1.)
88
89			! Calculer la longitude vraie de l'orbite terrestre:
90	guez	130	CALL orbite(real(jour), zlonsun, zdist)
91	guez	125
92			! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
93			zdeclin = asin(sin(zlonsunzpi/180.0)sin(r_incl*zpi/180.0))
94
95			DO i = 1, klon
96			aa = sin(rlat(i)zpi/180.0)sin(zdeclin)
97			bb = cos(rlat(i)zpi/180.0)cos(zdeclin)
98
99			! Midi local (angle du temps = 0.0):
100			rmu = aa + bb*cos(0.0)
101			rmu = max(0.0, rmu)
102			! IM cf. PB alb = 0.058/(rmu + 0.30)
103			! alb = 0.058/(rmu + 0.30) * 1.5
104			alb = 0.058/(rmu+0.30)*1.2
105			! alb = 0.058/(rmu + 0.30) * 1.3
106			srmu = rmu
107			salb = alb*rmu
108
109			! Faire l'integration numerique de midi a minuit (le facteur 2
110			! prend en compte l'autre moitie de la journee):
111			DO k = 1, npts
112			rmu = aa + bbcos(float(k)/float(npts)zpi)
113			rmu = max(0.0, rmu)
114			! IM cf. PB alb = 0.058/(rmu + 0.30)
115			! alb = 0.058/(rmu + 0.30) * 1.5
116			alb = 0.058/(rmu+0.30)*1.2
117			! alb = 0.058/(rmu + 0.30) * 1.3
118			srmu = srmu + rmu*2.0
119			salb = salb + albrmu2.0
120			END DO
121			IF (srmu/=0.0) THEN
122			albedo(i) = salb/srmu*fmagic
123			ELSE ! nuit polaire (on peut prendre une valeur quelconque)
124			albedo(i) = fmagic
125			END IF
126			END DO
127			END IF
128
129			END SUBROUTINE alboc
130
131			end module alboc_m