| 1 | SUBROUTINE aerave_new ( ndata, |
|---|
| 2 | & longdata,epdata,omegdata,gdata, |
|---|
| 3 | & longref,epref,temp,nir,longir |
|---|
| 4 | & ,epir,omegir,gir,qref,omegaref ) |
|---|
| 5 | |
|---|
| 6 | |
|---|
| 7 | IMPLICIT NONE |
|---|
| 8 | c....................................................................... |
|---|
| 9 | c |
|---|
| 10 | c R.Fournier 02/1996 |
|---|
| 11 | c (modif F.Forget 02/1996) |
|---|
| 12 | c le spectre est decoupe en "nir" bandes et cette routine calcule |
|---|
| 13 | c les donnees radiatives moyenne sur chaque bande : l'optimisation |
|---|
| 14 | c est faite pour une temperature au sol "temp" et une epaisseur |
|---|
| 15 | c optique de l'atmosphere "epref" a la longueur d'onde "longref" |
|---|
| 16 | c |
|---|
| 17 | c dans la version actuelle, les ponderations sont independantes de |
|---|
| 18 | c l'epaisseur optique : c'est a dire que "omegir", "gir" |
|---|
| 19 | c et "epir/epre" sont independants de "epref". |
|---|
| 20 | c en effet les ponderations sont choisies pour une solution exacte |
|---|
| 21 | c en couche mince et milieu isotherme. |
|---|
| 22 | c |
|---|
| 23 | c entree |
|---|
| 24 | c |
|---|
| 25 | c ndata : taille des champs data |
|---|
| 26 | c longdata,epdata,omegdata,gdata : proprietes radiative de l'aerosol |
|---|
| 27 | c (longdata longueur d'onde en METRES) |
|---|
| 28 | c * longref : longueur d'onde a laquelle l'epaisseur optique |
|---|
| 29 | c est connue |
|---|
| 30 | c * epref : epaisseur optique a longref |
|---|
| 31 | c * temp : temperature choisie pour la ponderation (Planck) |
|---|
| 32 | c * nir : nombre d'intervals dans la discretisation spectrale |
|---|
| 33 | c du GCM |
|---|
| 34 | c * longir : longueurs d'onde definissant ces intervals |
|---|
| 35 | c |
|---|
| 36 | c sortie |
|---|
| 37 | c |
|---|
| 38 | c * epir : epaisseur optique moyenne pour chaque interval |
|---|
| 39 | c * omegir : "scattering albedo" moyen pour chaque interval |
|---|
| 40 | c * gir : "assymetry factor" moyen pour chaque interval |
|---|
| 41 | c * qref : extinction coefficient at reference wavelength |
|---|
| 42 | c * omegaref : single scat. albedo at reference wavelength |
|---|
| 43 | c |
|---|
| 44 | c....................................................................... |
|---|
| 45 | c |
|---|
| 46 | REAL longref |
|---|
| 47 | REAL epref |
|---|
| 48 | REAL temp |
|---|
| 49 | INTEGER nir |
|---|
| 50 | REAL*8 longir(nir+1) |
|---|
| 51 | REAL epir(nir) |
|---|
| 52 | REAL omegir(nir) |
|---|
| 53 | REAL gir(nir) |
|---|
| 54 | c |
|---|
| 55 | c....................................................................... |
|---|
| 56 | c |
|---|
| 57 | INTEGER iir |
|---|
| 58 | INTEGER,PARAMETER :: nirmx=1900 |
|---|
| 59 | INTEGER idata,ndata |
|---|
| 60 | c |
|---|
| 61 | c....................................................................... |
|---|
| 62 | c |
|---|
| 63 | REAL emit |
|---|
| 64 | REAL totalemit(nirmx) |
|---|
| 65 | REAL longdata(ndata),epdata(ndata) |
|---|
| 66 | & ,omegdata(ndata),gdata(ndata) |
|---|
| 67 | REAL qextcorrdata(ndata) |
|---|
| 68 | INTEGER ibande |
|---|
| 69 | INTEGER,PARAMETER :: nbande=1000 |
|---|
| 70 | REAL long,deltalong |
|---|
| 71 | INTEGER ilong |
|---|
| 72 | INTEGER i1,i2 |
|---|
| 73 | REAL c1,c2 |
|---|
| 74 | REAL factep,qextcorr,omeg,g |
|---|
| 75 | REAL qref,omegaref |
|---|
| 76 | c |
|---|
| 77 | c....................................................................... |
|---|
| 78 | c |
|---|
| 79 | DOUBLE PRECISION tmp1 |
|---|
| 80 | REAL tmp2,tmp3 |
|---|
| 81 | c |
|---|
| 82 | c |
|---|
| 83 | long=longref |
|---|
| 84 | |
|---|
| 85 | c check ordering of longdata |
|---|
| 86 | DO idata=1,ndata-1 |
|---|
| 87 | IF (longdata(1).LT.longdata(ndata)) THEN |
|---|
| 88 | IF (.not.(longdata(idata).LT.longdata(idata+1))) THEN |
|---|
| 89 | call abort_physic("aerave_new", |
|---|
| 90 | & "Non descending order in longdata",1) |
|---|
| 91 | ENDIF |
|---|
| 92 | ELSEIF (longdata(1).GT.longdata(ndata)) THEN |
|---|
| 93 | IF (.not.(longdata(idata).GT.longdata(idata+1))) THEN |
|---|
| 94 | call abort_physic("aerave_new", |
|---|
| 95 | & "Non ascending order in longdata",1) |
|---|
| 96 | ENDIF |
|---|
| 97 | ENDIF |
|---|
| 98 | ENDDO |
|---|
| 99 | c |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | |
|---|
| 104 | c******************************************************** |
|---|
| 105 | c interpolation |
|---|
| 106 | c wavelengths (longdata) from data file in ascending order |
|---|
| 107 | IF (longdata(1).LT.longdata(ndata)) THEN |
|---|
| 108 | ilong=1 |
|---|
| 109 | DO idata=2,ndata |
|---|
| 110 | IF (long.gt.longdata(idata)) ilong=idata |
|---|
| 111 | ENDDO |
|---|
| 112 | i1=ilong |
|---|
| 113 | i2=ilong+1 |
|---|
| 114 | IF (i2.gt.ndata) i2=ndata |
|---|
| 115 | IF (long.lt.longdata(1)) i2=1 |
|---|
| 116 | IF (i1.eq.i2) THEN |
|---|
| 117 | c1=1.E+0 |
|---|
| 118 | c2=0.E+0 |
|---|
| 119 | ELSE |
|---|
| 120 | c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) |
|---|
| 121 | c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) |
|---|
| 122 | ENDIF |
|---|
| 123 | qref=c1*epdata(i1)+c2*epdata(i2) |
|---|
| 124 | omegaref=c1*omegdata(i1)+c2*omegdata(i2) |
|---|
| 125 | factep=qref/epref |
|---|
| 126 | DO idata=1,ndata |
|---|
| 127 | qextcorrdata(idata)=epdata(idata)/factep |
|---|
| 128 | ENDDO |
|---|
| 129 | c wavelengths (longdata) from data file in descending order |
|---|
| 130 | ELSEIF (longdata(1).GT.longdata(ndata)) THEN |
|---|
| 131 | ilong=1 |
|---|
| 132 | DO idata=2,ndata |
|---|
| 133 | IF (long.lt.longdata(idata)) ilong=idata |
|---|
| 134 | ENDDO |
|---|
| 135 | i1=ilong+1 |
|---|
| 136 | i2=ilong |
|---|
| 137 | IF (i1.gt.ndata) i1=ndata |
|---|
| 138 | IF (long.gt.longdata(1)) i1=1 |
|---|
| 139 | IF (i1.eq.i2) THEN |
|---|
| 140 | c1=1.E+0 |
|---|
| 141 | c2=0.E+0 |
|---|
| 142 | ELSE |
|---|
| 143 | c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) |
|---|
| 144 | c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) |
|---|
| 145 | ENDIF |
|---|
| 146 | qref=c1*epdata(i1)+c2*epdata(i2) |
|---|
| 147 | omegaref=c1*omegdata(i1)+c2*omegdata(i2) |
|---|
| 148 | factep=qref/epref |
|---|
| 149 | DO idata=1,ndata |
|---|
| 150 | qextcorrdata(idata)=epdata(idata)/factep |
|---|
| 151 | ENDDO |
|---|
| 152 | ENDIF |
|---|
| 153 | |
|---|
| 154 | c******************************************************** |
|---|
| 155 | c....................................................................... |
|---|
| 156 | c wavelengths (longdata) from data file in ascending order |
|---|
| 157 | c....................................................................... |
|---|
| 158 | IF (longdata(1).LT.longdata(ndata)) THEN |
|---|
| 159 | DO iir=1,nir |
|---|
| 160 | c |
|---|
| 161 | c....................................................................... |
|---|
| 162 | c |
|---|
| 163 | deltalong=(longir(iir+1)-longir(iir)) / nbande |
|---|
| 164 | totalemit(iir)=0.E+0 |
|---|
| 165 | epir(iir)=0.E+0 |
|---|
| 166 | omegir(iir)=0.E+0 |
|---|
| 167 | gir(iir)=0.E+0 |
|---|
| 168 | c |
|---|
| 169 | c....................................................................... |
|---|
| 170 | c |
|---|
| 171 | DO ibande=1,nbande |
|---|
| 172 | c |
|---|
| 173 | c....................................................................... |
|---|
| 174 | c |
|---|
| 175 | long=longir(iir) + (ibande-0.5E+0) * deltalong |
|---|
| 176 | CALL blackl(DBLE(long),DBLE(temp),tmp1) |
|---|
| 177 | emit=REAL(tmp1) |
|---|
| 178 | c |
|---|
| 179 | c....................................................................... |
|---|
| 180 | c |
|---|
| 181 | c interpolation |
|---|
| 182 | ilong=1 |
|---|
| 183 | DO idata=2,ndata |
|---|
| 184 | IF (long.gt.longdata(idata)) ilong=idata |
|---|
| 185 | ENDDO |
|---|
| 186 | i1=ilong |
|---|
| 187 | i2=ilong+1 |
|---|
| 188 | IF (i2.gt.ndata) i2=ndata |
|---|
| 189 | IF (long.lt.longdata(1)) i2=1 |
|---|
| 190 | IF (i1.eq.i2) THEN |
|---|
| 191 | c1=1.E+0 |
|---|
| 192 | c2=0.E+0 |
|---|
| 193 | ELSE |
|---|
| 194 | c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) |
|---|
| 195 | c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) |
|---|
| 196 | ENDIF |
|---|
| 197 | qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2) |
|---|
| 198 | omeg=c1*omegdata(i1)+c2*omegdata(i2) |
|---|
| 199 | g=c1*gdata(i1)+c2*gdata(i2) |
|---|
| 200 | c |
|---|
| 201 | c....................................................................... |
|---|
| 202 | c |
|---|
| 203 | totalemit(iir)=totalemit(iir)+deltalong*emit |
|---|
| 204 | epir(iir)=epir(iir)+deltalong*emit*qextcorr |
|---|
| 205 | omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr |
|---|
| 206 | gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g |
|---|
| 207 | c |
|---|
| 208 | c....................................................................... |
|---|
| 209 | c |
|---|
| 210 | ENDDO |
|---|
| 211 | c |
|---|
| 212 | c....................................................................... |
|---|
| 213 | c |
|---|
| 214 | gir(iir)=gir(iir)/omegir(iir) |
|---|
| 215 | omegir(iir)=omegir(iir)/epir(iir) |
|---|
| 216 | epir(iir)=epir(iir)/totalemit(iir) |
|---|
| 217 | c |
|---|
| 218 | c....................................................................... |
|---|
| 219 | c |
|---|
| 220 | ENDDO |
|---|
| 221 | c....................................................................... |
|---|
| 222 | c wavelengths (longdata) from data file in descending order |
|---|
| 223 | c....................................................................... |
|---|
| 224 | ELSEIF (longdata(1).GT.longdata(ndata)) THEN |
|---|
| 225 | DO iir=1,nir |
|---|
| 226 | c |
|---|
| 227 | c....................................................................... |
|---|
| 228 | c |
|---|
| 229 | deltalong=(longir(iir+1)-longir(iir)) / nbande |
|---|
| 230 | totalemit(iir)=0.E+0 |
|---|
| 231 | epir(iir)=0.E+0 |
|---|
| 232 | omegir(iir)=0.E+0 |
|---|
| 233 | gir(iir)=0.E+0 |
|---|
| 234 | c |
|---|
| 235 | c....................................................................... |
|---|
| 236 | c |
|---|
| 237 | DO ibande=1,nbande |
|---|
| 238 | c |
|---|
| 239 | c....................................................................... |
|---|
| 240 | c |
|---|
| 241 | long=longir(iir) + (ibande-0.5E+0) * deltalong |
|---|
| 242 | CALL blackl(DBLE(long),DBLE(temp),tmp1) |
|---|
| 243 | emit=REAL(tmp1) |
|---|
| 244 | c |
|---|
| 245 | c....................................................................... |
|---|
| 246 | c |
|---|
| 247 | c interpolation |
|---|
| 248 | ilong=1 |
|---|
| 249 | DO idata=2,ndata |
|---|
| 250 | IF (long.lt.longdata(idata)) ilong=idata |
|---|
| 251 | ENDDO |
|---|
| 252 | i1=ilong+1 |
|---|
| 253 | i2=ilong |
|---|
| 254 | IF (i1.gt.ndata) i1=ndata |
|---|
| 255 | IF (long.gt.longdata(1)) i1=1 |
|---|
| 256 | IF (i1.eq.i2) THEN |
|---|
| 257 | c1=1.E+0 |
|---|
| 258 | c2=0.E+0 |
|---|
| 259 | ELSE |
|---|
| 260 | c1=(longdata(i2)-long) / (longdata(i2)-longdata(i1)) |
|---|
| 261 | c2=(longdata(i1)-long) / (longdata(i1)-longdata(i2)) |
|---|
| 262 | ENDIF |
|---|
| 263 | qextcorr=c1*qextcorrdata(i1)+c2*qextcorrdata(i2) |
|---|
| 264 | omeg=c1*omegdata(i1)+c2*omegdata(i2) |
|---|
| 265 | g=c1*gdata(i1)+c2*gdata(i2) |
|---|
| 266 | c |
|---|
| 267 | c....................................................................... |
|---|
| 268 | c |
|---|
| 269 | totalemit(iir)=totalemit(iir)+deltalong*emit |
|---|
| 270 | epir(iir)=epir(iir)+deltalong*emit*qextcorr |
|---|
| 271 | omegir(iir)=omegir(iir)+deltalong*emit*omeg*qextcorr |
|---|
| 272 | gir(iir)=gir(iir)+deltalong*emit*omeg*qextcorr*g |
|---|
| 273 | c |
|---|
| 274 | c....................................................................... |
|---|
| 275 | c |
|---|
| 276 | ENDDO |
|---|
| 277 | c |
|---|
| 278 | c....................................................................... |
|---|
| 279 | c |
|---|
| 280 | gir(iir)=gir(iir)/omegir(iir) |
|---|
| 281 | omegir(iir)=omegir(iir)/epir(iir) |
|---|
| 282 | epir(iir)=epir(iir)/totalemit(iir) |
|---|
| 283 | c |
|---|
| 284 | c....................................................................... |
|---|
| 285 | c |
|---|
| 286 | ENDDO |
|---|
| 287 | ENDIF |
|---|
| 288 | c |
|---|
| 289 | c******************************************************** |
|---|
| 290 | c |
|---|
| 291 | c...................................................................... |
|---|
| 292 | c |
|---|
| 293 | c Diagnostic de controle si on moyenne sur tout le spectre vis ou IR : |
|---|
| 294 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 295 | c tmp2=0.E+0 |
|---|
| 296 | c DO iir=1,nir |
|---|
| 297 | c tmp2=tmp2+totalemit(iir) |
|---|
| 298 | c ENDDO |
|---|
| 299 | c tmp3=5.67E-8 * temp**4 |
|---|
| 300 | c IF (abs((tmp2-tmp3)/tmp3).gt.0.05E+0) THEN |
|---|
| 301 | c PRINT *,'!!!! <---> il manque du Planck (voir moyenne.F)' |
|---|
| 302 | c PRINT *,'somme des bandes :',tmp2,'--- Planck:',tmp3 |
|---|
| 303 | c ENDIF |
|---|
| 304 | c |
|---|
| 305 | c...................................................................... |
|---|
| 306 | c |
|---|
| 307 | END |
|---|
