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

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

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

trunk/libf/phylmd/albedo.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/Sources/phylmd/Interface_surf/alboc.f revision 208 by guez, Wed Dec 7 16:44:53 2016 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    pure function alboc(jour, rlat)
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      ! Author: Z. X. Li (LMD/CNRS)
12  c======================================================================      ! Date : 16 mars 1995
13  c Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)      ! Objet : Calculer l'alb\'edo sur l'oc\'ean
14  c Date: le 16 mars 1995      ! M\'ethode: int\'egrer num\'eriquement l'alb\'edo pendant une journ\'ee
15  c Objet: Calculer l'albedo sur l'ocean  
16  c Methode: Integrer numeriquement l'albedo pendant une journee      use nr_util, only: pi
17  c      USE orbite_m, ONLY: orbite
18  c Arguments;      USE yomcst, only: r_incl
19  c rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)  
20  c rlat (in,R)   : latitude en degre      integer, intent(in):: jour ! jour dans l'annee (a compter du 1 janvier)
21  c albedo (out,R): albedo obtenu (de 0 a 1)      REAL, intent(in):: rlat(:) ! latitude en degre
22  c======================================================================      real alboc(size(rlat)) ! albedo obtenu (de 0 a 1)
23  c  
24        REAL fmagic ! un facteur magique pour regler l'albedo      ! Local:
25  ccc      PARAMETER (fmagic=0.7)  
26  cccIM => a remplacer        INTEGER, PARAMETER:: npts =120
27  c       PARAMETER (fmagic=1.32)      ! Contr\^ole la pr\'ecision de l'int\'egration. 120 correspond \`a
28          PARAMETER (fmagic=1.0)      ! l'intervalle 6 minutes.
29  c       PARAMETER (fmagic=0.7)  
30        INTEGER npts ! il controle la precision de l'integration      REAL lonsun, declin
31        PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes      REAL rmu, srmu, salb, aa, bb
32  c      INTEGER i, k
33        REAL rlat(klon), rjour, albedo(klon)  
34        REAL zdist, zlonsun, zpi, zdeclin      !----------------------------------------------------------------------
35        REAL rmu,alb, srmu, salb, fauxo, aa, bb  
36        INTEGER i, k      ! Calculer la longitude vraie de l'orbite terrestre:
37  cccIM      CALL orbite(real(jour), lonsun)
38        LOGICAL ancien_albedo  
39        PARAMETER(ancien_albedo=.FALSE.)      ! Calculer la declinaison du soleil (qui varie entre + R_incl et - R_incl) :
40  c     SAVE albedo      declin = asin(sin(lonsun * pi / 180.) * sin(r_incl * pi / 180.))
41  c  
42        IF ( ancien_albedo ) THEN      DO i = 1, size(rlat)
43  c         aa = sin(rlat(i) * pi / 180.) * sin(declin)
44        zpi = 4. * ATAN(1.)         bb = cos(rlat(i) * pi / 180.) * cos(declin)
45  c  
46  c Calculer la longitude vraie de l'orbite terrestre:         ! Midi local (angle du temps = 0.):
47        CALL orbite(rjour,zlonsun,zdist)         rmu = max(0., aa + bb)
48  c         srmu = rmu
49  c Calculer la declinaison du soleil (qui varie entre + et - R_incl):         salb = 0.058 / (rmu + 0.30) * 1.2 * rmu
50        zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))  
51  c         ! Faire l'integration numerique de midi a minuit (le facteur 2
52        DO 999 i=1,klon         ! prend en compte l'autre moitie de la journee):
53        aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)         DO k = 1, npts
54        bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)            rmu = max(0., aa + bb * cos(real(k) / real(npts) * pi))
55  c            srmu = srmu + rmu * 2.
56  c Midi local (angle du temps = 0.0):            salb = salb + 0.058 / (rmu + 0.30) * 1.2 * rmu * 2.
57        rmu = aa + bb * COS(0.0)         END DO
58        rmu = MAX(0.0, rmu)         IF (srmu /= 0.) THEN
59        fauxo = (1.47-ACOS(rmu))/.15            alboc(i) = salb / srmu
60        alb = 0.03+0.630/(1.+fauxo*fauxo)         ELSE
61        srmu = rmu            ! nuit polaire (on peut prendre une valeur quelconque)
62        salb = alb * rmu            alboc(i) = 1.
63  c         END IF
64  c Faire l'integration numerique de midi a minuit (le facteur 2      END DO
65  c prend en compte l'autre moitie de la journee):  
66        DO k = 1, npts    END function alboc
67           rmu = aa + bb * COS(FLOAT(k)/FLOAT(npts)*zpi)  
68           rmu = MAX(0.0, rmu)  end module alboc_m
          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  
       ENDDO  
       IF (srmu .NE. 0.0) THEN  
          albedo(i) = salb / srmu * fmagic  
       ELSE ! nuit polaire (on peut prendre une valeur quelconque)  
          albedo(i) = fmagic  
       ENDIF  
   999 CONTINUE  
 c  
 c nouvel albedo  
 c  
       ELSE  
 c  
       zpi = 4. * ATAN(1.)  
 c  
 c Calculer la longitude vraie de l'orbite terrestre:  
       CALL orbite(rjour,zlonsun,zdist)  
 c  
 c Calculer la declinaison du soleil (qui varie entre + et - R_incl):  
       zdeclin = ASIN(SIN(zlonsun*zpi/180.0)*SIN(R_incl*zpi/180.0))  
 c  
       DO 1999 i=1,klon  
       aa = SIN(rlat(i)*zpi/180.0) * SIN(zdeclin)  
       bb = COS(rlat(i)*zpi/180.0) * COS(zdeclin)  
 c  
 c Midi local (angle du temps = 0.0):  
       rmu = aa + bb * COS(0.0)  
       rmu = MAX(0.0, rmu)  
 cIM cf. PB  alb = 0.058/(rmu + 0.30)  
 c     alb = 0.058/(rmu + 0.30) * 1.5  
       alb = 0.058/(rmu + 0.30) * 1.2  
 c     alb = 0.058/(rmu + 0.30) * 1.3  
       srmu = rmu  
       salb = alb * rmu  
 c  
 c Faire l'integration numerique de midi a minuit (le facteur 2  
 c 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)  
 cIM cf. PB      alb = 0.058/(rmu + 0.30)  
 c        alb = 0.058/(rmu + 0.30) * 1.5  
          alb = 0.058/(rmu + 0.30) * 1.2  
 c        alb = 0.058/(rmu + 0.30) * 1.3  
          srmu = srmu + rmu * 2.0  
          salb = salb + alb*rmu * 2.0  
       ENDDO  
       IF (srmu .NE. 0.0) THEN  
          albedo(i) = salb / srmu * fmagic  
       ELSE ! nuit polaire (on peut prendre une valeur quelconque)  
          albedo(i) = fmagic  
       ENDIF  
 1999  CONTINUE  
       ENDIF  
       RETURN  
       END  
 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.3  
changed lines
  Added in v.208
  ViewVC Help
Powered by ViewVC 1.1.21