/[lmdze]/trunk/phylmd/Interface_surf/alboc.f
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/alboc.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/phylmd/albedo.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/phylmd/Interface_surf/alboc.f revision 125 by guez, Fri Feb 6 15:00:28 2015 UTC
# Line 1  Line 1 
1  !  module alboc_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/albedo.F,v 1.2 2005/02/07 15:00:52 fairhead Exp $  
3  !    IMPLICIT NONE
4  c  
5  c  contains
6        SUBROUTINE alboc(rjour,rlat,albedo)  
7        use dimens_m    SUBROUTINE alboc(rjour, rlat, albedo)
8        use dimphy  
9        use yomcst      ! From LMDZ4/libf/phylmd/albedo.F, version 1.2 2005/02/07 15:00:52
10        use orbite_m, only: orbite  
11        IMPLICIT none      ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
12  c======================================================================      ! Date: le 16 mars 1995
13  c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)      ! Objet: Calculer l'albedo sur l'ocean
14  c Date: le 16 mars 1995      ! Methode: Integrer numeriquement l'albedo pendant une journee
15  c Objet: Calculer l'albedo sur l'ocean  
16  c Methode: Integrer numeriquement l'albedo pendant une journee      USE dimphy, only: klon
17  c      USE yomcst, only: r_incl
18  c Arguments;      USE orbite_m, ONLY: orbite
19  c rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)  
20  c rlat (in,R)   : latitude en degre      ! Arguments;
21  c albedo (out,R): albedo obtenu (de 0 a 1)      ! rjour (in, R) : jour dans l'annee (a compter du 1 janvier)
22  c======================================================================      ! rlat (in, R) : latitude en degre
23  c      ! albedo (out, R): albedo obtenu (de 0 a 1)
24        REAL fmagic ! un facteur magique pour regler l'albedo  
25  ccc      PARAMETER (fmagic=0.7)      REAL fmagic ! un facteur magique pour regler l'albedo
26  cccIM => a remplacer        ! cc PARAMETER (fmagic=0.7)
27  c       PARAMETER (fmagic=1.32)      ! ccIM => a remplacer
28          PARAMETER (fmagic=1.0)      ! PARAMETER (fmagic=1.32)
29  c       PARAMETER (fmagic=0.7)      PARAMETER (fmagic=1.0)
30        INTEGER npts ! il controle la precision de l'integration      ! PARAMETER (fmagic=0.7)
31        PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes      INTEGER npts ! il controle la precision de l'integration
32  c      PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes
33        REAL rlat(klon), rjour, albedo(klon)  
34        REAL zdist, zlonsun, zpi, zdeclin      REAL rlat(klon), rjour, albedo(klon)
35        REAL rmu,alb, srmu, salb, fauxo, aa, bb      REAL zdist, zlonsun, zpi, zdeclin
36        INTEGER i, k      REAL rmu, alb, srmu, salb, fauxo, aa, bb
37  cccIM      INTEGER i, k
38        LOGICAL ancien_albedo      ! ccIM
39        PARAMETER(ancien_albedo=.FALSE.)      LOGICAL ancien_albedo
40  c     SAVE albedo      PARAMETER (ancien_albedo=.FALSE.)
41  c      ! SAVE albedo
42        IF ( ancien_albedo ) THEN  
43  c      IF (ancien_albedo) THEN
44        zpi = 4. * ATAN(1.)  
45  c         zpi = 4.*atan(1.)
46  c Calculer la longitude vraie de l'orbite terrestre:  
47        CALL orbite(rjour,zlonsun,zdist)         ! Calculer la longitude vraie de l'orbite terrestre:
48  c         CALL orbite(rjour, zlonsun, zdist)
49  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):  
50        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))         ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
51  c         zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
52        DO 999 i=1,klon  
53        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)         DO i = 1, klon
54        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)            aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
55  c            bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
56  c Midi local (angle du temps = 0.0):  
57        rmu = aa + bb * COS(0.0)            ! Midi local (angle du temps = 0.0):
58        rmu = MAX(0.0, rmu)            rmu = aa + bb*cos(0.0)
59        fauxo = (1.47-ACOS(rmu))/.15            rmu = max(0.0, rmu)
60        alb = 0.03+0.630/(1.+fauxo*fauxo)            fauxo = (1.47-acos(rmu))/.15
61        srmu = rmu            alb = 0.03 + 0.630/(1.+fauxo*fauxo)
62        salb = alb * rmu            srmu = rmu
63  c            salb = alb*rmu
64  c Faire l'integration numerique de midi a minuit (le facteur 2  
65  c prend en compte l'autre moitie de la journee):            ! Faire l'integration numerique de midi a minuit (le facteur 2
66        DO k = 1, npts            ! prend en compte l'autre moitie de la journee):
67           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)            DO k = 1, npts
68           rmu = MAX(0.0, rmu)               rmu = aa + bb*cos(float(k)/float(npts)*zpi)
69           fauxo = (1.47-ACOS(rmu))/.15               rmu = max(0.0, rmu)
70           alb = 0.03+0.630/(1.+fauxo*fauxo)               fauxo = (1.47-acos(rmu))/.15
71           srmu = srmu + rmu * 2.0               alb = 0.03 + 0.630/(1.+fauxo*fauxo)
72           salb = salb + alb*rmu * 2.0               srmu = srmu + rmu*2.0
73        ENDDO               salb = salb + alb*rmu*2.0
74        IF (srmu .NE. 0.0) THEN            END DO
75           albedo(i) = salb / srmu * fmagic            IF (srmu/=0.0) THEN
76        ELSE ! nuit polaire (on peut prendre une valeur quelconque)               albedo(i) = salb/srmu*fmagic
77           albedo(i) = fmagic            ELSE ! nuit polaire (on peut prendre une valeur quelconque)
78        ENDIF               albedo(i) = fmagic
79    999 CONTINUE            END IF
80  c         END DO
81  c nouvel albedo  
82  c         ! nouvel albedo
83        ELSE  
84  c      ELSE
85        zpi = 4. * ATAN(1.)  
86  c         zpi = 4.*atan(1.)
87  c Calculer la longitude vraie de l'orbite terrestre:  
88        CALL orbite(rjour,zlonsun,zdist)         ! Calculer la longitude vraie de l'orbite terrestre:
89  c         CALL orbite(rjour, zlonsun, zdist)
90  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):  
91        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))         ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
92  c         zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
93        DO 1999 i=1,klon  
94        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)         DO i = 1, klon
95        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)            aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
96  c            bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
97  c Midi local (angle du temps = 0.0):  
98        rmu = aa + bb * COS(0.0)            ! Midi local (angle du temps = 0.0):
99        rmu = MAX(0.0, rmu)            rmu = aa + bb*cos(0.0)
100  cIM cf. PB  alb = 0.058/(rmu + 0.30)            rmu = max(0.0, rmu)
101  c     alb = 0.058/(rmu + 0.30) * 1.5            ! IM cf. PB alb = 0.058/(rmu + 0.30)
102        alb = 0.058/(rmu + 0.30) * 1.2            ! alb = 0.058/(rmu + 0.30) * 1.5
103  c     alb = 0.058/(rmu + 0.30) * 1.3            alb = 0.058/(rmu+0.30)*1.2
104        srmu = rmu            ! alb = 0.058/(rmu + 0.30) * 1.3
105        salb = alb * rmu            srmu = rmu
106  c            salb = alb*rmu
107  c Faire l'integration numerique de midi a minuit (le facteur 2  
108  c prend en compte l'autre moitie de la journee):            ! Faire l'integration numerique de midi a minuit (le facteur 2
109        DO k = 1, npts            ! prend en compte l'autre moitie de la journee):
110           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)            DO k = 1, npts
111           rmu = MAX(0.0, rmu)               rmu = aa + bb*cos(float(k)/float(npts)*zpi)
112  cIM cf. PB      alb = 0.058/(rmu + 0.30)               rmu = max(0.0, rmu)
113  c        alb = 0.058/(rmu + 0.30) * 1.5               ! IM cf. PB alb = 0.058/(rmu + 0.30)
114           alb = 0.058/(rmu + 0.30) * 1.2               ! alb = 0.058/(rmu + 0.30) * 1.5
115  c        alb = 0.058/(rmu + 0.30) * 1.3               alb = 0.058/(rmu+0.30)*1.2
116           srmu = srmu + rmu * 2.0               ! alb = 0.058/(rmu + 0.30) * 1.3
117           salb = salb + alb*rmu * 2.0               srmu = srmu + rmu*2.0
118        ENDDO               salb = salb + alb*rmu*2.0
119        IF (srmu .NE. 0.0) THEN            END DO
120           albedo(i) = salb / srmu * fmagic            IF (srmu/=0.0) THEN
121        ELSE ! nuit polaire (on peut prendre une valeur quelconque)               albedo(i) = salb/srmu*fmagic
122           albedo(i) = fmagic            ELSE ! nuit polaire (on peut prendre une valeur quelconque)
123        ENDIF               albedo(i) = fmagic
124  1999  CONTINUE            END IF
125        ENDIF         END DO
126        RETURN      END IF
127        END  
128  c=====================================================================    END SUBROUTINE alboc
129        SUBROUTINE alboc_cd(rmu0,albedo)  
130        use dimens_m  end module alboc_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========================================================================  
Legend:
Removed from v.76  
changed lines
  Added in v.125
  ViewVC Help
Powered by ViewVC 1.1.21