| 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=kflev) ! fichiers Vincent: same grid than GCM |
|---|
| 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 |
|---|
| 63 | real lambda(nnuve) ! wavelenght in table (mu->m, middle of interval) |
|---|
| 64 | real ksive(0:nlve+1,0:nlve+1,nnuve) ! ksi factors |
|---|
| 65 | real bplck(0:nlve+1,nnuve) ! Planck luminances in table layers |
|---|
| 66 | real al(nnuve),bl(nnuve) ! for Planck luminances calculations |
|---|
| 67 | real psive(0:nlve+1,0:nlve+1,nnuve) ! NER in W/m**2 per wavelength band |
|---|
| 68 | real psi_1(0:nlve+1,0:nlve+1) ! NER in W/m**2 (sum on lambda) |
|---|
| 69 | |
|---|
| 70 | real ztemp(0:nlve) ! GCM temperature in table layers |
|---|
| 71 | real zlnet(nlve+1) ! net thermal flux (W/m**2) |
|---|
| 72 | real dzlnet(0:nlve) ! Radiative budget (W/m**2) |
|---|
| 73 | character*22 nullchar |
|---|
| 74 | real lambdamin,lambdamax ! in microns |
|---|
| 75 | real dlambda ! cm-1 |
|---|
| 76 | |
|---|
| 77 | real y(0:nlve,nnuve) ! intermediaire Planck |
|---|
| 78 | real pdp(nlve) ! epaisseur de la couche en pression (Pa) |
|---|
| 79 | real zdblay(nlve,nnuve) ! gradient en temperature de planck |
|---|
| 80 | |
|---|
| 81 | logical firstcall |
|---|
| 82 | data firstcall/.true./ |
|---|
| 83 | |
|---|
| 84 | save lambda,ksive,al,bl,firstcall |
|---|
| 85 | |
|---|
| 86 | c ------------------------ |
|---|
| 87 | c Loading the files |
|---|
| 88 | c ------------------------ |
|---|
| 89 | |
|---|
| 90 | if (firstcall) then |
|---|
| 91 | |
|---|
| 92 | c Matrice Ksi |
|---|
| 93 | c------------ |
|---|
| 94 | open(13,file='ksi_gccr.txt') |
|---|
| 95 | read(13,*) nl0,nnu0 |
|---|
| 96 | if (nl0.ne.nlve) then |
|---|
| 97 | print*,'Probleme de dimension entre ksi.txt et lw' |
|---|
| 98 | print*,'N levels = ',nl0,nlve |
|---|
| 99 | stop |
|---|
| 100 | endif |
|---|
| 101 | if (nnu0.ne.nnuve) then |
|---|
| 102 | print*,'Probleme de dimension entre ksi.txt et lw' |
|---|
| 103 | print*,'N freq = ',nnu0,nnuve |
|---|
| 104 | stop |
|---|
| 105 | endif |
|---|
| 106 | do band=1,nnuve |
|---|
| 107 | read(13,*) lambdamin,lambdamax ! en microns |
|---|
| 108 | lambda(band)=(lambdamin+lambdamax)/2.*1.e-6 ! en m |
|---|
| 109 | dlambda =(1./lambdamin-1./lambdamax)*1.e4 ! en cm-1 |
|---|
| 110 | c print*,band,lambdamin,dlambda,lambdamax |
|---|
| 111 | do i=0,nlve+1 |
|---|
| 112 | read(13,'(52e17.9)') (ksive(i,j,band),j=0,nlve+1) |
|---|
| 113 | c ecart-type MC sur les ksi: pas utilise |
|---|
| 114 | c read(13,'(52e17.9)') (psive(i,j,band),j=0,nlve+1) |
|---|
| 115 | c changement de convention (signe) pour ksi, |
|---|
| 116 | c et prise en compte de la largeur de bande (en cm-1): |
|---|
| 117 | do j=0,nlve+1 |
|---|
| 118 | ksive(i,j,band) = -ksive(i,j,band)*dlambda |
|---|
| 119 | enddo |
|---|
| 120 | enddo |
|---|
| 121 | c calcul des coeff al et bl pour luminance Planck |
|---|
| 122 | al(band) = 2.*RHPLA*RCLUM*RCLUM/(lambda(band))**5. |
|---|
| 123 | c cette luminance doit etre en W/m²/sr/µm pour correspondre au calcul |
|---|
| 124 | c des ksi. Ici, elle est en W/m²/sr/m donc il faut mettre un facteur 1.e-6 |
|---|
| 125 | . * 1.e-6 |
|---|
| 126 | bl(band) = RHPLA*RCLUM/(RKBOL*lambda(band)) |
|---|
| 127 | enddo |
|---|
| 128 | close(13) |
|---|
| 129 | |
|---|
| 130 | endif ! firstcall |
|---|
| 131 | |
|---|
| 132 | c -------------------------------------- |
|---|
| 133 | c Calculation of the Psi matrix |
|---|
| 134 | c -------------------------------------- |
|---|
| 135 | |
|---|
| 136 | c temperature in the table layers |
|---|
| 137 | c ------------------------------- |
|---|
| 138 | |
|---|
| 139 | do j=1,nlve |
|---|
| 140 | ztemp(j) = PT(j) |
|---|
| 141 | enddo |
|---|
| 142 | |
|---|
| 143 | ztemp(0) = PTSURF |
|---|
| 144 | |
|---|
| 145 | c Planck function |
|---|
| 146 | c --------------- |
|---|
| 147 | |
|---|
| 148 | do band=1,nnuve |
|---|
| 149 | y(0,band) = exp(bl(band)/ztemp(0))-1. |
|---|
| 150 | bplck(0,band) = al(band)/(y(0,band)) |
|---|
| 151 | do j=1,nlve |
|---|
| 152 | c Developpement en polynomes ? |
|---|
| 153 | c bplck(j,band) = xp(1,band) |
|---|
| 154 | c . +ztemp(j)*(xp(2,band) |
|---|
| 155 | c . +ztemp(j)*(xp(3,band) |
|---|
| 156 | c . +ztemp(j)*(xp(4,band) |
|---|
| 157 | c . +ztemp(j)*(xp(5,band) |
|---|
| 158 | c . +ztemp(j)*(xp(6,band) ))))) |
|---|
| 159 | |
|---|
| 160 | c B(T,l) = al/(exp(bl/T)-1) |
|---|
| 161 | y(j,band) = exp(bl(band)/ztemp(j))-1. |
|---|
| 162 | bplck(j,band) = al(band)/(y(j,band)) |
|---|
| 163 | zdblay(j,band) = al(band)*bl(band)*exp(bl(band)/ztemp(j))/ |
|---|
| 164 | . ((ztemp(j)**2)*(y(j,band)**2)) |
|---|
| 165 | enddo |
|---|
| 166 | bplck(nlve+1,band) = 0.0 |
|---|
| 167 | enddo |
|---|
| 168 | |
|---|
| 169 | c Calculation of Psi |
|---|
| 170 | c ------------------ |
|---|
| 171 | |
|---|
| 172 | do band=1,nnuve |
|---|
| 173 | do j=0,nlve+1 |
|---|
| 174 | do i=0,nlve+1 |
|---|
| 175 | psive(i,j,band)=ksive(i,j,band)*(bplck(i,band)-bplck(j,band)) |
|---|
| 176 | enddo |
|---|
| 177 | enddo |
|---|
| 178 | enddo |
|---|
| 179 | |
|---|
| 180 | do j=0,nlve+1 |
|---|
| 181 | do i=0,nlve+1 |
|---|
| 182 | psi_1(i,j) = 0.0 ! positif quand nrj de i->j |
|---|
| 183 | enddo |
|---|
| 184 | enddo |
|---|
| 185 | |
|---|
| 186 | do band=1,nnuve |
|---|
| 187 | do j=0,nlve+1 |
|---|
| 188 | do i=0,nlve+1 |
|---|
| 189 | psi_1(i,j) = psi_1(i,j)+psive(i,j,band) |
|---|
| 190 | enddo |
|---|
| 191 | enddo |
|---|
| 192 | enddo |
|---|
| 193 | |
|---|
| 194 | c Verif...----------------------- |
|---|
| 195 | c open(11,file="psi.dat") |
|---|
| 196 | c do i=0,nlve+1 |
|---|
| 197 | c write(11,'(I3,83E17.9)') i,(psi_1(j,i),j=0,nlve+1) |
|---|
| 198 | c enddo |
|---|
| 199 | c close(11) |
|---|
| 200 | c stop |
|---|
| 201 | c ------------------------------- |
|---|
| 202 | |
|---|
| 203 | c -------------------------- |
|---|
| 204 | c Calculation of the fluxes |
|---|
| 205 | c -------------------------- |
|---|
| 206 | |
|---|
| 207 | c flux aux intercouches: |
|---|
| 208 | c zlnet(i+1) est le flux net traversant le plafond de la couche i (+ vers le haut) |
|---|
| 209 | do p=0,nlve ! numero de la couche |
|---|
| 210 | zlnet(p+1) = 0.0 |
|---|
| 211 | do j=p+1,nlve+1 |
|---|
| 212 | do i=0,p |
|---|
| 213 | zlnet(p+1) = zlnet(p+1)+psi_1(i,j) |
|---|
| 214 | enddo |
|---|
| 215 | enddo |
|---|
| 216 | enddo |
|---|
| 217 | |
|---|
| 218 | c flux net au sol, + vers le haut: |
|---|
| 219 | PSOLLW = zlnet(1) |
|---|
| 220 | c flux vers le bas au sol, + vers le bas: |
|---|
| 221 | PSOLLWDN = 0.0 |
|---|
| 222 | do i=1,nlve+1 |
|---|
| 223 | PSOLLWDN = PSOLLWDN+max(psi_1(i,0),0.0) |
|---|
| 224 | enddo |
|---|
| 225 | |
|---|
| 226 | c dfluxnet = radiative budget (W m-2) |
|---|
| 227 | do p=0,nlve ! numero de la couche |
|---|
| 228 | dzlnet(p) = 0.0 |
|---|
| 229 | do j=0,nlve+1 |
|---|
| 230 | dzlnet(p) = dzlnet(p)+psi_1(p,j) |
|---|
| 231 | enddo |
|---|
| 232 | enddo |
|---|
| 233 | |
|---|
| 234 | |
|---|
| 235 | c -------------------------------------- |
|---|
| 236 | c Interpolation in the GCM vertical grid |
|---|
| 237 | c -------------------------------------- |
|---|
| 238 | |
|---|
| 239 | c Flux net |
|---|
| 240 | c -------- |
|---|
| 241 | |
|---|
| 242 | do j=1,kflev+1 |
|---|
| 243 | ZFLNET(j) = zlnet(j) |
|---|
| 244 | enddo |
|---|
| 245 | |
|---|
| 246 | PTOPLW = ZFLNET(kflev+1) |
|---|
| 247 | |
|---|
| 248 | c Heating rates |
|---|
| 249 | c ------------- |
|---|
| 250 | |
|---|
| 251 | c cool (K/s) = dfluxnet (W/m2) ! positif quand nrj sort de la couche |
|---|
| 252 | c *g (m/s2) |
|---|
| 253 | c /(-dp) (epaisseur couche, en Pa=kg/m/s2) |
|---|
| 254 | c /cp (J/kg/K) |
|---|
| 255 | |
|---|
| 256 | |
|---|
| 257 | do j=1,kflev |
|---|
| 258 | pdp(j)=(PPB(j)-PPB(j+1))*1.e5 |
|---|
| 259 | enddo |
|---|
| 260 | |
|---|
| 261 | c calcul direct OU calcul par schema implicit |
|---|
| 262 | if (1.eq.0) then |
|---|
| 263 | do j=1,kflev |
|---|
| 264 | PCOOL(j) = dzlnet(j) *RG/RCPD / pdp(j) |
|---|
| 265 | enddo |
|---|
| 266 | else |
|---|
| 267 | call lwi(nlve,nnuve,dzlnet,zdblay,pdp,ksive,PCOOL) |
|---|
| 268 | endif |
|---|
| 269 | c print*,dzlnet |
|---|
| 270 | c print*,pdp |
|---|
| 271 | c print*,PCOOL |
|---|
| 272 | |
|---|
| 273 | do j=1,kflev |
|---|
| 274 | PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday |
|---|
| 275 | enddo |
|---|
| 276 | |
|---|
| 277 | firstcall = .false. |
|---|
| 278 | return |
|---|
| 279 | end |
|---|
| 280 | |
|---|
