| 1 | SUBROUTINE LW_venus_ve( |
|---|
| 2 | S PPB, PT, PTSURF, |
|---|
| 3 | S PCOOL, |
|---|
| 4 | S PTOPLW,PSOLLW,PSOLLWDN, |
|---|
| 5 | S ZFLNET) |
|---|
| 6 | |
|---|
| 7 | IMPLICIT none |
|---|
| 8 | |
|---|
| 9 | #include "dimensions.h" |
|---|
| 10 | #include "dimphy.h" |
|---|
| 11 | #include "raddim.h" |
|---|
| 12 | #include "YOMCST.h" |
|---|
| 13 | C |
|---|
| 14 | C ------------------------------------------------------------------ |
|---|
| 15 | C |
|---|
| 16 | C PURPOSE. |
|---|
| 17 | C -------- |
|---|
| 18 | C |
|---|
| 19 | c This routine loads the longwave matrix of factors Ksi, |
|---|
| 20 | c used to build the Net Exchange Rates matrix Psi. |
|---|
| 21 | c Psi(i,j,nu) = Ksi(i,j,nu) * ( B(i,nu)-B(j,nu) ) |
|---|
| 22 | c |
|---|
| 23 | c This Ksi matrix has been computed by Vincent Eymet |
|---|
| 24 | C |
|---|
| 25 | c The NER matrix is then integrated in frequency, and the output |
|---|
| 26 | c are calculated. |
|---|
| 27 | c |
|---|
| 28 | C AUTHOR. |
|---|
| 29 | C ------- |
|---|
| 30 | C Sebastien Lebonnois |
|---|
| 31 | C |
|---|
| 32 | C MODIFICATIONS. |
|---|
| 33 | C -------------- |
|---|
| 34 | C ORIGINAL : 27/07/2005 |
|---|
| 35 | C ------------------------------------------------------------------ |
|---|
| 36 | C |
|---|
| 37 | C* ARGUMENTS: |
|---|
| 38 | C |
|---|
| 39 | c inputs |
|---|
| 40 | |
|---|
| 41 | REAL PPB(KFLEV+1) ! inter-couches PRESSURE (bar) |
|---|
| 42 | REAL PT(KFLEV) ! Temperature in layer (K) |
|---|
| 43 | REAL PTSURF ! Surface temperature |
|---|
| 44 | C |
|---|
| 45 | c output |
|---|
| 46 | |
|---|
| 47 | REAL PCOOL(KFLEV) ! LONGWAVE COOLING (K/VENUSDAY) within each layer |
|---|
| 48 | REAL PTOPLW ! LONGWAVE FLUX AT T.O.A. (net, + vers le haut) |
|---|
| 49 | REAL PSOLLW ! LONGWAVE FLUX AT SURFACE (net, + vers le haut) |
|---|
| 50 | REAL PSOLLWDN ! LONGWAVE FLUX AT SURFACE (down, + vers le bas) |
|---|
| 51 | REAL ZFLNET(KFLEV+1) ! net thermal flux at ppb levels (+ vers le haut) |
|---|
| 52 | |
|---|
| 53 | C |
|---|
| 54 | C* LOCAL VARIABLES: |
|---|
| 55 | C |
|---|
| 56 | integer nlve,nnuve |
|---|
| 57 | parameter (nlve=81) ! fichiers Vincent |
|---|
| 58 | parameter (nnuve=68) ! fichiers Vincent et Bullock |
|---|
| 59 | real dureejour |
|---|
| 60 | parameter (dureejour=10.087e6) |
|---|
| 61 | |
|---|
| 62 | integer i,j,p,nl0,nnu0,band,k,l |
|---|
| 63 | real presve(nlve+1) ! pressure levels in table (Pa->bar) |
|---|
| 64 | real tempve(nlve+1) ! temperature in table (K) (middle of layer) |
|---|
| 65 | real altve(nlve+1) ! altitude in table (km) |
|---|
| 66 | real lambda(nnuve) ! wavelenght in table (mu->m, middle of interval) |
|---|
| 67 | real ksive(0:nlve+1,0:nlve+1,nnuve) ! ksi factors |
|---|
| 68 | real bplck(0:nlve+1,nnuve) ! Planck luminances in table layers |
|---|
| 69 | real al(nnuve),bl(nnuve) ! for Planck luminances calculations |
|---|
| 70 | real psive(0:nlve+1,0:nlve+1,nnuve) ! NER in W/m**2 per wavelength band |
|---|
| 71 | real psi_1(0:nlve+1,0:nlve+1) ! NER in W/m**2 (sum on lambda) |
|---|
| 72 | |
|---|
| 73 | real ztemp(0:nlve) ! GCM temperature in table layers |
|---|
| 74 | real zlnet(nlve+1) ! net thermal flux (W/m**2) |
|---|
| 75 | real dzlnet(0:nlve) ! Radiative budget (W/m**2) |
|---|
| 76 | real radbudget(kflev) ! Radiative budget on GCM grid |
|---|
| 77 | real coolrate(nlve) ! cooling rates (K/s) on table grid |
|---|
| 78 | character*22 nullchar |
|---|
| 79 | real lambdamin,lambdamax ! in microns |
|---|
| 80 | real dlambda ! cm-1 |
|---|
| 81 | |
|---|
| 82 | real y(0:nlve,nnuve) ! intermediaire Planck |
|---|
| 83 | real pdp(kflev) ! delta pression (Pa), grille GCM |
|---|
| 84 | real pdpve(nlve) ! delta pression (Pa), grille table |
|---|
| 85 | real zdblay(nlve,nnuve) ! gradient en temperature de planck |
|---|
| 86 | |
|---|
| 87 | real factflux |
|---|
| 88 | real facttemp,prT(kflev),prTve(nlve) |
|---|
| 89 | |
|---|
| 90 | logical firstcall |
|---|
| 91 | data firstcall/.true./ |
|---|
| 92 | |
|---|
| 93 | save lambda,ksive,al,bl,firstcall |
|---|
| 94 | |
|---|
| 95 | c ------------------------ |
|---|
| 96 | c Loading the files |
|---|
| 97 | c ------------------------ |
|---|
| 98 | |
|---|
| 99 | if (firstcall) then |
|---|
| 100 | |
|---|
| 101 | print*,"PREMIER APPEL RADIATIF" |
|---|
| 102 | |
|---|
| 103 | c Grilles alt et press |
|---|
| 104 | c--------------------- |
|---|
| 105 | open(11,file='mesh.txt') |
|---|
| 106 | read(11,*) nl0,nnu0,i |
|---|
| 107 | read(11,*) nullchar |
|---|
| 108 | read(11,'(82(2x,F15.9))') altve |
|---|
| 109 | read(11,'(82(2x,F15.9))') presve |
|---|
| 110 | read(11,'(81(2x,F15.9))') tempve |
|---|
| 111 | tempve(nlve+1)=tempve(nlve) |
|---|
| 112 | close(11) |
|---|
| 113 | if (nl0.ne.nlve) then |
|---|
| 114 | print*,'Probleme de dimension entre mesh.txt et lw' |
|---|
| 115 | print*,'N levels = ',nl0,nlve |
|---|
| 116 | stop |
|---|
| 117 | endif |
|---|
| 118 | if (nnu0.ne.nnuve) then |
|---|
| 119 | print*,'Probleme de dimension entre mesh.txt et lw' |
|---|
| 120 | print*,'N freq = ',nnu0,nnuve |
|---|
| 121 | stop |
|---|
| 122 | endif |
|---|
| 123 | do i=1,nlve+1 |
|---|
| 124 | presve(i) = presve(i)*1.e-5 ! convert to bar |
|---|
| 125 | enddo |
|---|
| 126 | |
|---|
| 127 | c Verifs... |
|---|
| 128 | c print*, altve |
|---|
| 129 | c print*, presve |
|---|
| 130 | |
|---|
| 131 | c Matrice Ksi |
|---|
| 132 | c------------ |
|---|
| 133 | open(13,file='ksi_gccr.txt') |
|---|
| 134 | read(13,*) nl0,nnu0 |
|---|
| 135 | if (nl0.ne.nlve) then |
|---|
| 136 | print*,'Probleme de dimension entre ksi.txt et lw' |
|---|
| 137 | print*,'N levels = ',nl0,nlve |
|---|
| 138 | stop |
|---|
| 139 | endif |
|---|
| 140 | if (nnu0.ne.nnuve) then |
|---|
| 141 | print*,'Probleme de dimension entre ksi.txt et lw' |
|---|
| 142 | print*,'N freq = ',nnu0,nnuve |
|---|
| 143 | stop |
|---|
| 144 | endif |
|---|
| 145 | do band=1,nnuve |
|---|
| 146 | read(13,*) lambdamin,lambdamax ! en microns |
|---|
| 147 | lambda(band)=(lambdamin+lambdamax)/2.*1.e-6 ! en m |
|---|
| 148 | dlambda =(1./lambdamin-1./lambdamax)*1.e4 ! en cm-1 |
|---|
| 149 | c print*,band,lambdamin,dlambda,lambdamax |
|---|
| 150 | do i=0,nlve+1 |
|---|
| 151 | read(13,'(83e17.9)') (ksive(i,j,band),j=0,nlve+1) |
|---|
| 152 | c ecart-type MC sur les ksi: pas utilise |
|---|
| 153 | c read(13,'(83e17.9)') (psive(i,j,band),j=0,nlve+1) |
|---|
| 154 | c changement de convention (signe) pour ksi, |
|---|
| 155 | c et prise en compte de la largeur de bande (en cm-1): |
|---|
| 156 | do j=0,nlve+1 |
|---|
| 157 | ksive(i,j,band) = -ksive(i,j,band)*dlambda |
|---|
| 158 | enddo |
|---|
| 159 | enddo |
|---|
| 160 | c calcul des coeff al et bl pour luminance Planck |
|---|
| 161 | al(band) = 2.*RHPLA*RCLUM*RCLUM/(lambda(band))**5. |
|---|
| 162 | c cette luminance doit etre en W/m²/sr/µm pour correspondre au calcul |
|---|
| 163 | c des ksi. Ici, elle est en W/m²/sr/m donc il faut mettre un facteur 1.e-6 |
|---|
| 164 | . * 1.e-6 |
|---|
| 165 | bl(band) = RHPLA*RCLUM/(RKBOL*lambda(band)) |
|---|
| 166 | enddo |
|---|
| 167 | close(13) |
|---|
| 168 | |
|---|
| 169 | endif ! firstcall |
|---|
| 170 | |
|---|
| 171 | c -------------------------------------- |
|---|
| 172 | c Calculation of the Psi matrix |
|---|
| 173 | c -------------------------------------- |
|---|
| 174 | |
|---|
| 175 | c temperature in the table layers |
|---|
| 176 | c ------------------------------- |
|---|
| 177 | |
|---|
| 178 | do i=1,kflev |
|---|
| 179 | prT(i) = (PPB(i)+PPB(i+1))/2. |
|---|
| 180 | c prT(i) = 10.**((log10(PPB(i))+log10(PPB(i+1)))/2.) |
|---|
| 181 | enddo |
|---|
| 182 | |
|---|
| 183 | do j=1,nlve |
|---|
| 184 | prTve(j) = (presve(j)+presve(j+1))/2. |
|---|
| 185 | c prTve(j) = max(10.**((log10(presve(j))+log10(presve(j+1)))/2.) |
|---|
| 186 | c . ,1.e-5) |
|---|
| 187 | enddo |
|---|
| 188 | |
|---|
| 189 | do j=1,nlve |
|---|
| 190 | nl0 = 2 |
|---|
| 191 | do i=1,kflev-1 |
|---|
| 192 | if (prT(i).ge.prTve(j)) then |
|---|
| 193 | nl0 = i+1 |
|---|
| 194 | endif |
|---|
| 195 | enddo |
|---|
| 196 | |
|---|
| 197 | facttemp = (log10(prTve(j))-log10(prT(nl0-1))) |
|---|
| 198 | . /(log10(prT(nl0))-log10(prT(nl0-1))) |
|---|
| 199 | ztemp(j) = facttemp *PT(nl0) |
|---|
| 200 | . + (1.-facttemp)*PT(nl0-1) |
|---|
| 201 | |
|---|
| 202 | c write(100,*) prTve(j),ztemp(j) |
|---|
| 203 | enddo |
|---|
| 204 | |
|---|
| 205 | ztemp(0) = PTSURF |
|---|
| 206 | |
|---|
| 207 | c do j=1,kflev |
|---|
| 208 | c write(101,*) prT(j),PT(j) |
|---|
| 209 | c enddo |
|---|
| 210 | |
|---|
| 211 | c print*,'VERIF TEMP' |
|---|
| 212 | c print*,PTSURF,PT |
|---|
| 213 | c print*,ztemp |
|---|
| 214 | c print*,tempve |
|---|
| 215 | |
|---|
| 216 | c Planck function |
|---|
| 217 | c --------------- |
|---|
| 218 | |
|---|
| 219 | do band=1,nnuve |
|---|
| 220 | y(0,band) = exp(bl(band)/ztemp(0))-1. |
|---|
| 221 | bplck(0,band) = al(band)/(y(0,band)) |
|---|
| 222 | do j=1,nlve |
|---|
| 223 | c Developpement en polynomes ? |
|---|
| 224 | c bplck(j,band) = xp(1,band) |
|---|
| 225 | c . +ztemp(j)*(xp(2,band) |
|---|
| 226 | c . +ztemp(j)*(xp(3,band) |
|---|
| 227 | c . +ztemp(j)*(xp(4,band) |
|---|
| 228 | c . +ztemp(j)*(xp(5,band) |
|---|
| 229 | c . +ztemp(j)*(xp(6,band) ))))) |
|---|
| 230 | |
|---|
| 231 | c B(T,l) = al/(exp(bl/T)-1) |
|---|
| 232 | y(j,band) = exp(bl(band)/ztemp(j))-1. |
|---|
| 233 | bplck(j,band) = al(band)/(y(j,band)) |
|---|
| 234 | zdblay(j,band) = al(band)*bl(band)*exp(bl(band)/ztemp(j))/ |
|---|
| 235 | . ((ztemp(j)**2)*(y(j,band)**2)) |
|---|
| 236 | enddo |
|---|
| 237 | bplck(nlve+1,band) = 0.0 |
|---|
| 238 | enddo |
|---|
| 239 | |
|---|
| 240 | c Calculation of Psi |
|---|
| 241 | c ------------------ |
|---|
| 242 | |
|---|
| 243 | do band=1,nnuve |
|---|
| 244 | do j=0,nlve+1 |
|---|
| 245 | do i=0,nlve+1 |
|---|
| 246 | psive(i,j,band)=ksive(i,j,band)*(bplck(i,band)-bplck(j,band)) |
|---|
| 247 | enddo |
|---|
| 248 | enddo |
|---|
| 249 | enddo |
|---|
| 250 | |
|---|
| 251 | do j=0,nlve+1 |
|---|
| 252 | do i=0,nlve+1 |
|---|
| 253 | psi_1(i,j) = 0.0 ! positif quand nrj de i->j |
|---|
| 254 | enddo |
|---|
| 255 | enddo |
|---|
| 256 | |
|---|
| 257 | do band=1,nnuve |
|---|
| 258 | do j=0,nlve+1 |
|---|
| 259 | do i=0,nlve+1 |
|---|
| 260 | psi_1(i,j) = psi_1(i,j)+psive(i,j,band) |
|---|
| 261 | enddo |
|---|
| 262 | enddo |
|---|
| 263 | enddo |
|---|
| 264 | |
|---|
| 265 | c Verif...----------------------- |
|---|
| 266 | c open(11,file="psi.dat") |
|---|
| 267 | c do i=0,nlve+1 |
|---|
| 268 | c write(11,'(I3,83E17.9)') i,(psi_1(j,i),j=0,nlve+1) |
|---|
| 269 | c enddo |
|---|
| 270 | c close(11) |
|---|
| 271 | c stop |
|---|
| 272 | c ------------------------------- |
|---|
| 273 | |
|---|
| 274 | c -------------------------- |
|---|
| 275 | c Calculation of the fluxes |
|---|
| 276 | c -------------------------- |
|---|
| 277 | |
|---|
| 278 | c flux aux intercouches: |
|---|
| 279 | c zlnet(i+1) est le flux net traversant le plafond de la couche i (+ vers le haut) |
|---|
| 280 | do p=0,nlve ! numero de la couche |
|---|
| 281 | zlnet(p+1) = 0.0 |
|---|
| 282 | do j=p+1,nlve+1 |
|---|
| 283 | do i=0,p |
|---|
| 284 | zlnet(p+1) = zlnet(p+1)+psi_1(i,j) |
|---|
| 285 | enddo |
|---|
| 286 | enddo |
|---|
| 287 | enddo |
|---|
| 288 | |
|---|
| 289 | c do p=1,nlve |
|---|
| 290 | c write(102,*) presve(p),zlnet(p), |
|---|
| 291 | c . (zlnet(p+1)-zlnet(p))/(presve(p)-presve(p+1)) |
|---|
| 292 | c enddo |
|---|
| 293 | |
|---|
| 294 | c flux net au sol, + vers le haut: |
|---|
| 295 | PSOLLW = zlnet(1) |
|---|
| 296 | c flux vers le bas au sol, + vers le bas: |
|---|
| 297 | PSOLLWDN = 0.0 |
|---|
| 298 | do i=1,nlve+1 |
|---|
| 299 | PSOLLWDN = PSOLLWDN+max(psi_1(i,0),0.0) |
|---|
| 300 | enddo |
|---|
| 301 | |
|---|
| 302 | c dfluxnet = radiative budget (W m-2) |
|---|
| 303 | do p=0,nlve ! numero de la couche |
|---|
| 304 | dzlnet(p) = 0.0 |
|---|
| 305 | do j=0,nlve+1 |
|---|
| 306 | dzlnet(p) = dzlnet(p)+psi_1(p,j) |
|---|
| 307 | enddo |
|---|
| 308 | enddo |
|---|
| 309 | |
|---|
| 310 | |
|---|
| 311 | c -------------------------------------- |
|---|
| 312 | c Interpolation in the GCM vertical grid |
|---|
| 313 | c -------------------------------------- |
|---|
| 314 | |
|---|
| 315 | c Flux net |
|---|
| 316 | c -------- |
|---|
| 317 | |
|---|
| 318 | do j=1,kflev+1 |
|---|
| 319 | nl0 = 2 |
|---|
| 320 | do i=1,nlve |
|---|
| 321 | if (presve(i).ge.PPB(j)) then |
|---|
| 322 | nl0 = i+1 |
|---|
| 323 | endif |
|---|
| 324 | enddo |
|---|
| 325 | |
|---|
| 326 | factflux = (log10(max(PPB(j),presve(nlve+1))) |
|---|
| 327 | . -log10(presve(nl0-1))) |
|---|
| 328 | . /(log10(presve(nl0))-log10(presve(nl0-1))) |
|---|
| 329 | ZFLNET(j) = factflux *zlnet(nl0) |
|---|
| 330 | . + (1.-factflux)*zlnet(nl0-1) |
|---|
| 331 | |
|---|
| 332 | enddo |
|---|
| 333 | |
|---|
| 334 | PTOPLW = ZFLNET(kflev+1) |
|---|
| 335 | |
|---|
| 336 | c Heating rates |
|---|
| 337 | c ------------- |
|---|
| 338 | |
|---|
| 339 | c cool (K/s) = dfluxnet (W/m2) ! positif quand nrj sort de la couche |
|---|
| 340 | c *g (m/s2) |
|---|
| 341 | c /(-dp) (epaisseur couche, en Pa=kg/m/s2) |
|---|
| 342 | c /cp (J/kg/K) |
|---|
| 343 | |
|---|
| 344 | c layers thickness on each pressure grid (in Pa) |
|---|
| 345 | |
|---|
| 346 | do j=1,kflev |
|---|
| 347 | pdp(j)=(PPB(j)-PPB(j+1))*1.e5 |
|---|
| 348 | enddo |
|---|
| 349 | |
|---|
| 350 | do j=1,nlve |
|---|
| 351 | pdpve(j)=(presve(j)-presve(j+1))*1.e5 |
|---|
| 352 | enddo |
|---|
| 353 | |
|---|
| 354 | c CHOIX CALCUL DIRECT OU IMPLICIT |
|---|
| 355 | |
|---|
| 356 | c Ici, le budget radiatif est en calcul direct. |
|---|
| 357 | c On ne fait rien. Si on veut l'implicit, on autorise le test suivant: |
|---|
| 358 | |
|---|
| 359 | if (1.eq.0) then |
|---|
| 360 | |
|---|
| 361 | c Pour calcul par schema implicite, on obtient en sortie de lwi le coolrate. |
|---|
| 362 | c Donc on actualise le dzlnet par dzlnet=coolrate*(cp/g)*pdpve |
|---|
| 363 | |
|---|
| 364 | call lwi(nlve,nnuve,dzlnet,zdblay,pdpve,ksive,coolrate) |
|---|
| 365 | do j=1,nlve |
|---|
| 366 | dzlnet(j) = coolrate(j) *RCPD/RG *pdpve(j) |
|---|
| 367 | enddo |
|---|
| 368 | |
|---|
| 369 | endif |
|---|
| 370 | |
|---|
| 371 | c Interpolation on GCM grid of radiative budgets (dzlnet) |
|---|
| 372 | |
|---|
| 373 | c on divise l'energie deposee dans la couche par l'epaisseur |
|---|
| 374 | c on moyenne ensuite ces valeurs (creneaux sur grille VE) |
|---|
| 375 | c entre les niveaux de la grille GCM, et on multiplie ensuite par |
|---|
| 376 | c l'epaisseur (nouvelle grille) pour avoir l'energie deposee dans les |
|---|
| 377 | c couches GCM. |
|---|
| 378 | |
|---|
| 379 | i=1 |
|---|
| 380 | do j=1,kflev |
|---|
| 381 | if (PPB(j+1).ge.presve(i+1)) then |
|---|
| 382 | radbudget(j) = dzlnet(i)/(log10(presve(i+1))-log10(presve(i))) |
|---|
| 383 | . *(log10(PPB(j+1))-log10(PPB(j))) |
|---|
| 384 | else |
|---|
| 385 | l=i+1 |
|---|
| 386 | do while ((PPB(j+1).lt.presve(l+1)).and.(l.ne.nlve)) |
|---|
| 387 | l=l+1 |
|---|
| 388 | enddo |
|---|
| 389 | radbudget(j) = dzlnet(i)/(log10(presve(i+1))-log10(presve(i)))* |
|---|
| 390 | . (log10(presve(i+1))-log10(PPB(j))) |
|---|
| 391 | . +dzlnet(l)/(log10(presve(l+1))-log10(presve(l)))* |
|---|
| 392 | . (log10(PPB(j+1))-log10(presve(l))) |
|---|
| 393 | do k=i+2,l |
|---|
| 394 | radbudget(j) = radbudget(j)+dzlnet(k-1) |
|---|
| 395 | enddo |
|---|
| 396 | i=l |
|---|
| 397 | endif |
|---|
| 398 | enddo |
|---|
| 399 | |
|---|
| 400 | c do i=1,kflev |
|---|
| 401 | c print*,radbudget(i),prT(i) |
|---|
| 402 | c enddo |
|---|
| 403 | c do i=1,nlve |
|---|
| 404 | c print*,dzlnet(i),prTve(i) |
|---|
| 405 | c enddo |
|---|
| 406 | c stop |
|---|
| 407 | |
|---|
| 408 | c On obtient le coolrate en calculant: PCOOL = radbudget*(g/cp)/pdp |
|---|
| 409 | |
|---|
| 410 | do j=1,kflev |
|---|
| 411 | PCOOL(j) = radbudget(j) *RG/RCPD / pdp(j) |
|---|
| 412 | PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday |
|---|
| 413 | enddo |
|---|
| 414 | |
|---|
| 415 | c print*,PCOOL |
|---|
| 416 | |
|---|
| 417 | firstcall = .false. |
|---|
| 418 | return |
|---|
| 419 | end |
|---|
