[3] | 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 | |
---|