Changeset 5144 for LMDZ6/branches/Amaury_dev/libf/phylmd/albedo.F90
- Timestamp:
- Jul 29, 2024, 11:01:04 PM (3 months ago)
- File:
-
- 1 edited
-
LMDZ6/branches/Amaury_dev/libf/phylmd/albedo.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/albedo.F90
r5137 r5144 9 9 SUBROUTINE alboc(rjour, rlat, albedo) 10 10 USE dimphy 11 USE lmdz_yomcst 12 11 13 ! ====================================================================== 12 14 ! Auteur(s): Z.X. Li (LMD/CNRS) (adaptation du GCM du LMD) … … 20 22 ! albedo (out,R): albedo obtenu (de 0 a 1) 21 23 ! ====================================================================== 22 include "YOMCST.h"23 24 24 25 INTEGER npts ! il controle la precision de l'integration 25 PARAMETER (npts =120) ! 120 correspond a l'interval 6 minutes26 PARAMETER (npts = 120) ! 120 correspond a l'interval 6 minutes 26 27 27 28 REAL rlat(klon), rjour, albedo(klon) … … 31 32 ! ccIM 32 33 LOGICAL ancien_albedo 33 PARAMETER (ancien_albedo =.FALSE.)34 PARAMETER (ancien_albedo = .FALSE.) 34 35 ! SAVE albedo 35 36 36 37 IF (ancien_albedo) THEN 37 38 38 zpi = 4.*atan(1.)39 zpi = 4. * atan(1.) 39 40 40 41 41 ! Calculer la longitude vraie de l'orbite terrestre: 42 CALL orbite(rjour, zlonsun, zdist) 42 43 43 44 zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))44 ! Calculer la declinaison du soleil (qui varie entre + et - R_incl): 45 zdeclin = asin(sin(zlonsun * zpi / 180.0) * sin(r_incl * zpi / 180.0)) 45 46 46 47 aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)48 bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)47 DO i = 1, klon 48 aa = sin(rlat(i) * zpi / 180.0) * sin(zdeclin) 49 bb = cos(rlat(i) * zpi / 180.0) * cos(zdeclin) 49 50 50 ! Midi local (angle du temps = 0.0): 51 rmu = aa + bb*cos(0.0) 51 ! Midi local (angle du temps = 0.0): 52 rmu = aa + bb * cos(0.0) 53 rmu = max(0.0, rmu) 54 fauxo = (1.47 - acos(rmu)) / .15 55 alb = 0.03 + 0.630 / (1. + fauxo * fauxo) 56 srmu = rmu 57 salb = alb * rmu 58 59 ! Faire l'integration numerique de midi a minuit (le facteur 2 60 ! prend en compte l'autre moitie de la journee): 61 DO k = 1, npts 62 rmu = aa + bb * cos(real(k) / real(npts) * zpi) 52 63 rmu = max(0.0, rmu) 53 fauxo = (1.47-acos(rmu))/.15 54 alb = 0.03 + 0.630/(1.+fauxo*fauxo) 55 srmu = rmu 56 salb = alb*rmu 64 fauxo = (1.47 - acos(rmu)) / .15 65 alb = 0.03 + 0.630 / (1. + fauxo * fauxo) 66 srmu = srmu + rmu * 2.0 67 salb = salb + alb * rmu * 2.0 68 END DO 69 IF (srmu/=0.0) THEN 70 albedo(i) = salb / srmu 71 ELSE ! nuit polaire (on peut prendre une valeur quelconque) 72 albedo(i) = 1.0 73 END IF 74 END DO 57 75 58 ! Faire l'integration numerique de midi a minuit (le facteur 2 59 ! prend en compte l'autre moitie de la journee): 60 DO k = 1, npts 61 rmu = aa + bb*cos(real(k)/real(npts)*zpi) 62 rmu = max(0.0, rmu) 63 fauxo = (1.47-acos(rmu))/.15 64 alb = 0.03 + 0.630/(1.+fauxo*fauxo) 65 srmu = srmu + rmu*2.0 66 salb = salb + alb*rmu*2.0 67 END DO 68 IF (srmu/=0.0) THEN 69 albedo(i) = salb/srmu 70 ELSE ! nuit polaire (on peut prendre une valeur quelconque) 71 albedo(i) = 1.0 72 END IF 73 END DO 74 75 ! nouvel albedo 76 ! nouvel albedo 76 77 77 78 ELSE 78 79 79 zpi = 4.*atan(1.)80 zpi = 4. * atan(1.) 80 81 81 82 82 ! Calculer la longitude vraie de l'orbite terrestre: 83 CALL orbite(rjour, zlonsun, zdist) 83 84 84 85 zdeclin = asin(sin(zlonsun*zpi/180.0)*sin(r_incl*zpi/180.0))85 ! Calculer la declinaison du soleil (qui varie entre + et - R_incl): 86 zdeclin = asin(sin(zlonsun * zpi / 180.0) * sin(r_incl * zpi / 180.0)) 86 87 87 88 aa = sin(rlat(i)*zpi/180.0)*sin(zdeclin)89 bb = cos(rlat(i)*zpi/180.0)*cos(zdeclin)88 DO i = 1, klon 89 aa = sin(rlat(i) * zpi / 180.0) * sin(zdeclin) 90 bb = cos(rlat(i) * zpi / 180.0) * cos(zdeclin) 90 91 91 ! Midi local (angle du temps = 0.0): 92 rmu = aa + bb*cos(0.0) 92 ! Midi local (angle du temps = 0.0): 93 rmu = aa + bb * cos(0.0) 94 rmu = max(0.0, rmu) 95 ! IM cf. PB alb = 0.058/(rmu + 0.30) 96 ! alb = 0.058/(rmu + 0.30) * 1.5 97 alb = 0.058 / (rmu + 0.30) * 1.2 98 ! alb = 0.058/(rmu + 0.30) * 1.3 99 srmu = rmu 100 salb = alb * rmu 101 102 ! Faire l'integration numerique de midi a minuit (le facteur 2 103 ! prend en compte l'autre moitie de la journee): 104 DO k = 1, npts 105 rmu = aa + bb * cos(real(k) / real(npts) * zpi) 93 106 rmu = max(0.0, rmu) 94 ! IM cf. PB alb = 0.058/(rmu + 0.30)107 ! IM cf. PB alb = 0.058/(rmu + 0.30) 95 108 ! alb = 0.058/(rmu + 0.30) * 1.5 96 alb = 0.058 /(rmu+0.30)*1.2109 alb = 0.058 / (rmu + 0.30) * 1.2 97 110 ! alb = 0.058/(rmu + 0.30) * 1.3 98 srmu = rmu 99 salb = alb*rmu 100 101 ! Faire l'integration numerique de midi a minuit (le facteur 2 102 ! prend en compte l'autre moitie de la journee): 103 DO k = 1, npts 104 rmu = aa + bb*cos(real(k)/real(npts)*zpi) 105 rmu = max(0.0, rmu) 106 ! IM cf. PB alb = 0.058/(rmu + 0.30) 107 ! alb = 0.058/(rmu + 0.30) * 1.5 108 alb = 0.058/(rmu+0.30)*1.2 109 ! alb = 0.058/(rmu + 0.30) * 1.3 110 srmu = srmu + rmu*2.0 111 salb = salb + alb*rmu*2.0 112 END DO 113 IF (srmu/=0.0) THEN 114 albedo(i) = salb/srmu 115 ELSE ! nuit polaire (on peut prendre une valeur quelconque) 116 albedo(i) = 1.0 117 END IF 118 END DO 111 srmu = srmu + rmu * 2.0 112 salb = salb + alb * rmu * 2.0 113 END DO 114 IF (srmu/=0.0) THEN 115 albedo(i) = salb / srmu 116 ELSE ! nuit polaire (on peut prendre une valeur quelconque) 117 albedo(i) = 1.0 118 END IF 119 END DO 119 120 END IF 120 121 … … 136 137 ! albedo (OUT): albedo de surface de l'ocean 137 138 ! ====================================================================== 138 REAL, INTENT(IN) :: rmu0(klon)139 REAL, INTENT(OUT) :: albedo(klon)139 REAL, INTENT(IN) :: rmu0(klon) 140 REAL, INTENT(OUT) :: albedo(klon) 140 141 141 142 REAL fauxo 142 143 INTEGER i 143 144 LOGICAL ancien_albedo 144 PARAMETER (ancien_albedo =.FALSE.)145 PARAMETER (ancien_albedo = .FALSE.) 145 146 146 147 IF (ancien_albedo) THEN 147 148 fauxo = (1.47-acos(max(rmu0(i), 0.0)))/0.15149 albedo(i) = 0.03+.630/(1.+fauxo*fauxo)150 albedo(i) = max(min(albedo(i),0.60), 0.04)151 148 DO i = 1, klon 149 fauxo = (1.47 - acos(max(rmu0(i), 0.0))) / 0.15 150 albedo(i) = 0.03 + .630 / (1. + fauxo * fauxo) 151 albedo(i) = max(min(albedo(i), 0.60), 0.04) 152 END DO 152 153 ELSE 153 154 albedo(i) = 0.058/(max(rmu0(i), 0.0)+0.30)155 albedo(i) = max(min(albedo(i),0.60), 0.04)156 154 DO i = 1, klon 155 albedo(i) = 0.058 / (max(rmu0(i), 0.0) + 0.30) 156 albedo(i) = max(min(albedo(i), 0.60), 0.04) 157 END DO 157 158 END IF 158 159
Note: See TracChangeset
for help on using the changeset viewer.
