phylmd/Interface_surf/alboc.f

module alboc_m

  IMPLICIT NONE

contains

  pure function alboc(jour, rlat)

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

    ! Author: Z. X. Li (LMD/CNRS)
    ! Date : 16 mars 1995
    ! Objet : Calculer l'alb\'edo sur l'oc\'ean
    ! M\'ethode: int\'egrer num\'eriquement l'alb\'edo pendant une journ\'ee

    use nr_util, only: pi
    USE orbite_m, ONLY: orbite
    USE yomcst, only: r_incl

    integer, intent(in):: jour ! jour dans l'annee (a compter du 1 janvier)
    REAL, intent(in):: rlat(:) ! latitude en degre
    real alboc(size(rlat)) ! albedo obtenu (de 0 a 1)

    ! Local:

    INTEGER, PARAMETER:: npts =120
    ! Contr\^ole la pr\'ecision de l'int\'egration. 120 correspond \`a
    ! l'intervalle 6 minutes.

    REAL lonsun, declin
    REAL rmu, srmu, salb, aa, bb
    INTEGER i, k

    !----------------------------------------------------------------------

    ! Calculer la longitude vraie de l'orbite terrestre:
    CALL orbite(real(jour), lonsun)

    ! Calculer la declinaison du soleil (qui varie entre + R_incl et - R_incl) :
    declin = asin(sin(lonsun * pi / 180.) * sin(r_incl * pi / 180.))

    DO i = 1, size(rlat)
       aa = sin(rlat(i) * pi / 180.) * sin(declin)
       bb = cos(rlat(i) * pi / 180.) * cos(declin)

       ! Midi local (angle du temps = 0.):
       rmu = max(0., aa + bb)
       srmu = rmu
       salb = 0.058 / (rmu + 0.30) * 1.2 * 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 = max(0., aa + bb * cos(real(k) / real(npts) * pi))
          srmu = srmu + rmu * 2.
          salb = salb + 0.058 / (rmu + 0.30) * 1.2 * rmu * 2.
       END DO
       IF (srmu /= 0.) THEN
          alboc(i) = salb / srmu
       ELSE
          ! nuit polaire (on peut prendre une valeur quelconque)
          alboc(i) = 1.
       END IF
    END DO

  END function 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	208	pure function alboc(jour, rlat)
8	guez	81
9	guez	208	! From LMDZ4/libf/phylmd/albedo.F, version 1.2, 2005/02/07 15:00:52
10	guez	81
11	guez	208	! Author: Z. X. Li (LMD/CNRS)
12	guez	154	! Date : 16 mars 1995
13			! Objet : Calculer l'alb\'edo sur l'oc\'ean
14			! M\'ethode: int\'egrer num\'eriquement l'alb\'edo pendant une journ\'ee
15	guez	81
16	guez	154	use nr_util, only: pi
17			USE orbite_m, ONLY: orbite
18	guez	125	USE yomcst, only: r_incl
19	guez	81
20	guez	154	integer, intent(in):: jour ! jour dans l'annee (a compter du 1 janvier)
21			REAL, intent(in):: rlat(:) ! latitude en degre
22	guez	208	real alboc(size(rlat)) ! albedo obtenu (de 0 a 1)
23	guez	81
24	guez	154	! Local:
25	guez	81
26	guez	154	INTEGER, PARAMETER:: npts =120
27			! Contr\^ole la pr\'ecision de l'int\'egration. 120 correspond \`a
28			! l'intervalle 6 minutes.
29	guez	81
30	guez	208	REAL lonsun, declin
31			REAL rmu, srmu, salb, aa, bb
32	guez	154	INTEGER i, k
33	guez	81
34	guez	154	!----------------------------------------------------------------------
35	guez	81
36	guez	154	! Calculer la longitude vraie de l'orbite terrestre:
37	guez	208	CALL orbite(real(jour), lonsun)
38	guez	81
39	guez	208	! Calculer la declinaison du soleil (qui varie entre + R_incl et - R_incl) :
40			declin = asin(sin(lonsun * pi / 180.) * sin(r_incl * pi / 180.))
41	guez	81
42	guez	154	DO i = 1, size(rlat)
43	guez	208	aa = sin(rlat(i) * pi / 180.) * sin(declin)
44			bb = cos(rlat(i) * pi / 180.) * cos(declin)
45	guez	81
46	guez	208	! Midi local (angle du temps = 0.):
47			rmu = max(0., aa + bb)
48	guez	154	srmu = rmu
49	guez	208	salb = 0.058 / (rmu + 0.30) * 1.2 * rmu
50	guez	81
51	guez	154	! Faire l'integration numerique de midi a minuit (le facteur 2
52			! prend en compte l'autre moitie de la journee):
53			DO k = 1, npts
54	guez	208	rmu = max(0., aa + bb * cos(real(k) / real(npts) * pi))
55			srmu = srmu + rmu * 2.
56			salb = salb + 0.058 / (rmu + 0.30) * 1.2 * rmu * 2.
57	guez	125	END DO
58	guez	208	IF (srmu /= 0.) THEN
59			alboc(i) = salb / srmu
60			ELSE
61			! nuit polaire (on peut prendre une valeur quelconque)
62			alboc(i) = 1.
63	guez	154	END IF
64			END DO
65	guez	125
66	guez	208	END function alboc
67	guez	125
68			end module alboc_m