Changeset 2322


Ignore:
Timestamp:
Jul 7, 2015, 6:01:01 PM (9 years ago)
Author:
lguez
Message:

Bug fix in alboc_cd (more than 11 years old). rmu0 must not be
modified in alboc_cd. The corresponding actual argument in subroutine
surf_ocean is an intent(in) argument of surf_ocean. It probably worked
before because the compiler made a copy of the argument but I would
say the behavior was in principle compiler-dependent and
optimization-level-dependent.

Location:
LMDZ5/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/albedo.F90

    r1992 r2322  
     1! $Id$
     2module albedo
    13
    2 ! $Id$
     4  IMPLICIT NONE
    35
     6contains
    47
     8  SUBROUTINE alboc(rjour, rlat, albedo)
     9    USE dimphy
     10    ! ======================================================================
     11    ! 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
    515
    6 SUBROUTINE alboc(rjour, rlat, albedo)
    7   USE dimphy
    8   IMPLICIT NONE
    9   ! ======================================================================
    10   ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD)
    11   ! Date: le 16 mars 1995
    12   ! Objet: Calculer l'albedo sur l'ocean
    13   ! Methode: Integrer numeriquement l'albedo pendant une journee
     16    ! Arguments;
     17    ! rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)
     18    ! rlat (in,R)   : latitude en degre
     19    ! albedo (out,R): albedo obtenu (de 0 a 1)
     20    ! ======================================================================
     21    ! ym#include "dimensions.h"
     22    ! ym#include "dimphy.h"
     23    include "YOMCST.h"
     24    include "clesphys.h"
    1425
    15   ! Arguments;
    16   ! rjour (in,R)  : jour dans l'annee (a compter du 1 janvier)
    17   ! rlat (in,R)   : latitude en degre
    18   ! albedo (out,R): albedo obtenu (de 0 a 1)
    19   ! ======================================================================
    20   ! ym#include "dimensions.h"
    21   ! ym#include "dimphy.h"
    22   include "YOMCST.h"
    23   include "clesphys.h"
     26    ! fmagic -> clesphys.h/.inc
     27    ! REAL fmagic ! un facteur magique pour regler l'albedo
     28    ! cc      PARAMETER (fmagic=0.7)
     29    ! ccIM => a remplacer
     30    ! PARAMETER (fmagic=1.32)
     31    ! PARAMETER (fmagic=1.0)
     32    ! PARAMETER (fmagic=0.7)
     33    INTEGER npts ! il controle la precision de l'integration
     34    PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes
    2435
    25   ! fmagic -> clesphys.h/.inc
    26   ! REAL fmagic ! un facteur magique pour regler l'albedo
    27   ! cc      PARAMETER (fmagic=0.7)
    28   ! ccIM => a remplacer
    29   ! PARAMETER (fmagic=1.32)
    30   ! PARAMETER (fmagic=1.0)
    31   ! PARAMETER (fmagic=0.7)
    32   INTEGER npts ! il controle la precision de l'integration
    33   PARAMETER (npts=120) ! 120 correspond a l'interval 6 minutes
     36    REAL rlat(klon), rjour, albedo(klon)
     37    REAL zdist, zlonsun, zpi, zdeclin
     38    REAL rmu, alb, srmu, salb, fauxo, aa, bb
     39    INTEGER i, k
     40    ! ccIM
     41    LOGICAL ancien_albedo
     42    PARAMETER (ancien_albedo=.FALSE.)
     43    ! SAVE albedo
    3444
    35   REAL rlat(klon), rjour, albedo(klon)
    36   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
     45    IF (ancien_albedo) THEN
    4346
    44   IF (ancien_albedo) THEN
     47       zpi = 4.*atan(1.)
    4548
    46     zpi = 4.*atan(1.)
     49       ! Calculer la longitude vraie de l'orbite terrestre:
     50       CALL orbite(rjour, zlonsun, zdist)
    4751
    48     ! Calculer la longitude vraie de l'orbite terrestre:
    49     CALL orbite(rjour, zlonsun, zdist)
     52       ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
     53       zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
    5054
    51     ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
    52     zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
     55       DO i = 1, klon
     56          aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
     57          bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
    5358
    54     DO i = 1, klon
    55       aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
    56       bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
     59          ! Midi local (angle du temps = 0.0):
     60          rmu = aa + bb*cos(0.0)
     61          rmu = max(0.0, rmu)
     62          fauxo = (1.47-acos(rmu))/.15
     63          alb = 0.03 + 0.630/(1.+fauxo*fauxo)
     64          srmu = rmu
     65          salb = alb*rmu
    5766
    58       ! 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
     67          ! Faire l'integration numerique de midi a minuit (le facteur 2
     68          ! prend en compte l'autre moitie de la journee):
     69          DO k = 1, npts
     70             rmu = aa + bb*cos(real(k)/real(npts)*zpi)
     71             rmu = max(0.0, rmu)
     72             fauxo = (1.47-acos(rmu))/.15
     73             alb = 0.03 + 0.630/(1.+fauxo*fauxo)
     74             srmu = srmu + rmu*2.0
     75             salb = salb + alb*rmu*2.0
     76          END DO
     77          IF (srmu/=0.0) THEN
     78             albedo(i) = salb/srmu*fmagic + pmagic
     79          ELSE ! nuit polaire (on peut prendre une valeur quelconque)
     80             albedo(i) = fmagic
     81          END IF
     82       END DO
    6583
    66       ! 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 + bb*cos(real(k)/real(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 + alb*rmu*2.0
    75       END DO
    76       IF (srmu/=0.0) THEN
    77         albedo(i) = salb/srmu*fmagic + pmagic
    78       ELSE ! nuit polaire (on peut prendre une valeur quelconque)
    79         albedo(i) = fmagic
    80       END IF
    81     END DO
     84       ! nouvel albedo
    8285
    83     ! nouvel albedo
     86    ELSE
    8487
    85   ELSE
     88       zpi = 4.*atan(1.)
    8689
    87     zpi = 4.*atan(1.)
     90       ! Calculer la longitude vraie de l'orbite terrestre:
     91       CALL orbite(rjour, zlonsun, zdist)
    8892
    89     ! Calculer la longitude vraie de l'orbite terrestre:
    90     CALL orbite(rjour, zlonsun, zdist)
     93       ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
     94       zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
    9195
    92     ! Calculer la declinaison du soleil (qui varie entre + et - R_incl):
    93     zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))
     96       DO i = 1, klon
     97          aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)
     98          bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)
    9499
    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)
     100          ! Midi local (angle du temps = 0.0):
     101          rmu = aa + bb*cos(0.0)
     102          rmu = max(0.0, rmu)
     103          ! IM cf. PB  alb = 0.058/(rmu + 0.30)
     104          ! alb = 0.058/(rmu + 0.30) * 1.5
     105          alb = 0.058/(rmu+0.30)*1.2
     106          ! alb = 0.058/(rmu + 0.30) * 1.3
     107          srmu = rmu
     108          salb = alb*rmu
    98109
    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
     110          ! Faire l'integration numerique de midi a minuit (le facteur 2
     111          ! prend en compte l'autre moitie de la journee):
     112          DO k = 1, npts
     113             rmu = aa + bb*cos(real(k)/real(npts)*zpi)
     114             rmu = max(0.0, rmu)
     115             ! IM cf. PB      alb = 0.058/(rmu + 0.30)
     116             ! alb = 0.058/(rmu + 0.30) * 1.5
     117             alb = 0.058/(rmu+0.30)*1.2
     118             ! alb = 0.058/(rmu + 0.30) * 1.3
     119             srmu = srmu + rmu*2.0
     120             salb = salb + alb*rmu*2.0
     121          END DO
     122          IF (srmu/=0.0) THEN
     123             albedo(i) = salb/srmu*fmagic + pmagic
     124          ELSE ! nuit polaire (on peut prendre une valeur quelconque)
     125             albedo(i) = fmagic
     126          END IF
     127       END DO
     128    END IF
     129    RETURN
     130  END SUBROUTINE alboc
     131  ! =====================================================================
     132  SUBROUTINE alboc_cd(rmu0, albedo)
     133    USE dimphy
    108134
    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 + bb*cos(real(k)/real(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 + alb*rmu*2.0
    120       END DO
    121       IF (srmu/=0.0) THEN
    122         albedo(i) = salb/srmu*fmagic + pmagic
    123       ELSE ! nuit polaire (on peut prendre une valeur quelconque)
    124         albedo(i) = fmagic
    125       END IF
    126     END DO
    127   END IF
    128   RETURN
    129 END SUBROUTINE alboc
    130 ! =====================================================================
    131 SUBROUTINE alboc_cd(rmu0, albedo)
    132   USE dimphy
    133   IMPLICIT NONE
    134   ! ======================================================================
    135   ! Auteur(s): Z.X. Li (LMD/CNRS)
    136   ! date: 19940624
    137   ! Calculer l'albedo sur l'ocean en fonction de l'angle zenithal moyen
    138   ! Formule due a Larson and Barkstrom (1977) Proc. of the symposium
    139   ! on radiation in the atmosphere, 19-28 August 1976, science Press,
    140   ! 1977 pp 451-453, ou These de 3eme cycle de Sylvie Joussaume.
     135    ! Auteur(s): Z.X. Li (LMD/CNRS)
     136    ! date: 19940624
     137    ! Calculer l'albedo sur l'ocean en fonction de l'angle zenithal moyen
     138    ! Formule due a Larson and Barkstrom (1977) Proc. of the symposium
     139    ! on radiation in the atmosphere, 19-28 August 1976, science Press,
     140    ! 1977 pp 451-453, ou These de 3eme cycle de Sylvie Joussaume.
    141141
    142   ! Arguments
    143   ! rmu0    (in): cosinus de l'angle solaire zenithal
    144   ! albedo (out): albedo de surface de l'ocean
    145   ! ======================================================================
    146   ! ym#include "dimensions.h"
    147   ! ym#include "dimphy.h"
    148   include "clesphys.h"
    149   REAL rmu0(klon), albedo(klon)
     142    ! Arguments
     143    ! rmu0    (in): cosinus de l'angle solaire zenithal
     144    ! albedo (out): albedo de surface de l'ocean
     145    ! ======================================================================
     146    include "clesphys.h"
     147    REAL, intent(in):: rmu0(klon)
     148    real, intent(out):: albedo(klon)
    150149
    151   ! REAL fmagic ! un facteur magique pour regler l'albedo
    152   ! cc      PARAMETER (fmagic=0.7)
    153   ! ccIM => a remplacer
    154   ! PARAMETER (fmagic=1.32)
    155   ! PARAMETER (fmagic=1.0)
    156   ! PARAMETER (fmagic=0.7)
     150    ! REAL fmagic ! un facteur magique pour regler l'albedo
     151    ! cc      PARAMETER (fmagic=0.7)
     152    ! ccIM => a remplacer
     153    ! PARAMETER (fmagic=1.32)
     154    ! PARAMETER (fmagic=1.0)
     155    ! PARAMETER (fmagic=0.7)
    157156
    158   REAL fauxo
    159   INTEGER i
    160   ! ccIM
    161   LOGICAL ancien_albedo
    162   PARAMETER (ancien_albedo=.FALSE.)
    163   ! SAVE albedo
     157    REAL fauxo
     158    INTEGER i
     159    LOGICAL ancien_albedo
     160    PARAMETER (ancien_albedo=.FALSE.)
    164161
    165   IF (ancien_albedo) THEN
     162    IF (ancien_albedo) THEN
     163       DO i = 1, klon
     164          fauxo = (1.47-acos(max(rmu0(i), 0.0)))/0.15
     165          albedo(i) = fmagic*(.03+.630/(1.+fauxo*fauxo)) + pmagic
     166          albedo(i) = max(min(albedo(i),0.60), 0.04)
     167       END DO
     168    ELSE
     169       DO i = 1, klon
     170          albedo(i) = fmagic*0.058/(max(rmu0(i), 0.0)+0.30) + pmagic
     171          albedo(i) = max(min(albedo(i),0.60), 0.04)
     172       END DO
     173    END IF
    166174
    167     DO i = 1, klon
     175  END SUBROUTINE alboc_cd
    168176
    169       rmu0(i) = max(rmu0(i), 0.0)
    170 
    171       fauxo = (1.47-acos(rmu0(i)))/0.15
    172       albedo(i) = fmagic*(.03+.630/(1.+fauxo*fauxo)) + pmagic
    173       albedo(i) = max(min(albedo(i),0.60), 0.04)
    174     END DO
    175 
    176     ! nouvel albedo
    177 
    178   ELSE
    179 
    180     DO i = 1, klon
    181       rmu0(i) = max(rmu0(i), 0.0)
    182       ! IM:orig albedo(i) = 0.058/(rmu0(i) + 0.30)
    183       albedo(i) = fmagic*0.058/(rmu0(i)+0.30) + pmagic
    184       albedo(i) = max(min(albedo(i),0.60), 0.04)
    185     END DO
    186 
    187   END IF
    188 
    189   RETURN
    190 END SUBROUTINE alboc_cd
    191 ! ========================================================================
     177end module albedo

  • LMDZ5/trunk/libf/phylmd/surf_ocean_mod.F90

    r2278 r2322  
    2020       flux_u1, flux_v1)
    2121
     22  use albedo, only: alboc, alboc_cd
    2223  USE dimphy
    2324  USE surface_data, ONLY     : type_ocean
Note: See TracChangeset for help on using the changeset viewer.