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