| 1 | ! |
|---|
| 2 | ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $ |
|---|
| 3 | ! |
|---|
| 4 | SUBROUTINE radlwsw(dist, rmu0, fract, zzlev, |
|---|
| 5 | . paprs, pplay,tsol, t, |
|---|
| 6 | . heat,cool,radsol, |
|---|
| 7 | . topsw,toplw,solsw,sollw, |
|---|
| 8 | . sollwdown, |
|---|
| 9 | . lwnet, swnet) |
|---|
| 10 | c |
|---|
| 11 | c====================================================================== |
|---|
| 12 | c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 |
|---|
| 13 | c Objet: interface entre le modele et les rayonnements |
|---|
| 14 | c Arguments: |
|---|
| 15 | c dist-----input-R- distance astronomique terre-soleil |
|---|
| 16 | c rmu0-----input-R- cosinus de l'angle zenithal |
|---|
| 17 | c fract----input-R- duree d'ensoleillement normalisee |
|---|
| 18 | c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h) |
|---|
| 19 | c zzlev----input-R- altitude a inter-couche (m) |
|---|
| 20 | c paprs----input-R- pression a inter-couche (Pa) |
|---|
| 21 | c pplay----input-R- pression au milieu de couche (Pa) |
|---|
| 22 | c tsol-----input-R- temperature du sol (en K) |
|---|
| 23 | c t--------input-R- temperature (K) |
|---|
| 24 | c heat-----output-R- echauffement atmospherique (visible) (K/s) |
|---|
| 25 | c cool-----output-R- refroidissement dans l'IR (K/s) |
|---|
| 26 | c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas) |
|---|
| 27 | c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas) |
|---|
| 28 | c toplw----output-R- ray. IR net au sommet de l'atmosphere (+ vers le haut) |
|---|
| 29 | c solsw----output-R- flux solaire net a la surface (+ vers le bas) |
|---|
| 30 | c sollw----output-R- ray. IR net a la surface (+ vers le bas) |
|---|
| 31 | c sollwdown-output-R- ray. IR descendant a la surface (+ vers le bas) |
|---|
| 32 | c lwnet____output-R- flux IR net (+ vers le haut) |
|---|
| 33 | c swnet____output-R- flux solaire net (+ vers le bas) |
|---|
| 34 | c |
|---|
| 35 | |
|---|
| 36 | c MODIFS pour multimatrices ksi SPECIFIQUE VENUS |
|---|
| 37 | c S. Lebonnois 20/12/2006 |
|---|
| 38 | c corrections 13/07/2007 |
|---|
| 39 | c New ksi matrix: possibility of different cloud model fct of lat 05/2014 |
|---|
| 40 | |
|---|
| 41 | c====================================================================== |
|---|
| 42 | use dimphy |
|---|
| 43 | USE comgeomphy |
|---|
| 44 | use write_field_phy |
|---|
| 45 | IMPLICIT none |
|---|
| 46 | #include "dimensions.h" |
|---|
| 47 | #include "YOMCST.h" |
|---|
| 48 | #include "clesphys.h" |
|---|
| 49 | #include "comcstVE.h" |
|---|
| 50 | |
|---|
| 51 | !=========== |
|---|
| 52 | ! Arguments |
|---|
| 53 | !=========== |
|---|
| 54 | real rmu0(klon), fract(klon), dist |
|---|
| 55 | |
|---|
| 56 | REAL zzlev(klon,klev+1) |
|---|
| 57 | real paprs(klon,klev+1), pplay(klon,klev) |
|---|
| 58 | real tsol(klon) |
|---|
| 59 | real t(klon,klev) |
|---|
| 60 | real heat(klon,klev), cool(klon,klev) |
|---|
| 61 | real radsol(klon), topsw(klon), toplw(klon) |
|---|
| 62 | real solsw(klon), sollw(klon) |
|---|
| 63 | real sollwdown(klon) |
|---|
| 64 | REAL swnet(klon,klev+1),lwnet(klon,klev+1) |
|---|
| 65 | |
|---|
| 66 | !=========== |
|---|
| 67 | ! Local |
|---|
| 68 | !=========== |
|---|
| 69 | INTEGER k, kk, i, j, band |
|---|
| 70 | |
|---|
| 71 | REAL PPB(klev+1) |
|---|
| 72 | |
|---|
| 73 | REAL zfract, zrmu0 |
|---|
| 74 | |
|---|
| 75 | REAL zheat(klev), zcool(klev) |
|---|
| 76 | real temp(klev),znivs(klev+1) |
|---|
| 77 | REAL ZFSNET(klev+1),ZFLNET(klev+1) |
|---|
| 78 | REAL ztopsw, ztoplw |
|---|
| 79 | REAL zsolsw, zsollw |
|---|
| 80 | cIM BEG |
|---|
| 81 | REAL zsollwdown |
|---|
| 82 | cIM END |
|---|
| 83 | real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file |
|---|
| 84 | |
|---|
| 85 | real psi(0:klev+1,0:klev+1) |
|---|
| 86 | real deltapsi(0:klev+1,0:klev+1) |
|---|
| 87 | real pt0(0:klev+1) |
|---|
| 88 | real bplck(0:klev+1,nnuve) ! Planck luminances in table layers |
|---|
| 89 | real y(0:klev,nnuve) ! temporary variable for Planck |
|---|
| 90 | real zdblay(0:klev+1,nnuve) ! temperature gradient of planck function |
|---|
| 91 | integer mat0,lat,ips,isza,ips0,isza0 |
|---|
| 92 | real factp,factz,ksi |
|---|
| 93 | |
|---|
| 94 | logical firstcall |
|---|
| 95 | data firstcall/.true./ |
|---|
| 96 | save firstcall |
|---|
| 97 | |
|---|
| 98 | c------------------------------------------- |
|---|
| 99 | c Initialisations |
|---|
| 100 | c----------------- |
|---|
| 101 | |
|---|
| 102 | if (firstcall) then |
|---|
| 103 | |
|---|
| 104 | c ---------- ksive -------------- |
|---|
| 105 | allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat)) |
|---|
| 106 | call load_ksi(ksive) |
|---|
| 107 | |
|---|
| 108 | endif ! firstcall |
|---|
| 109 | c------------------------------------------- |
|---|
| 110 | |
|---|
| 111 | DO k = 1, klev |
|---|
| 112 | DO i = 1, klon |
|---|
| 113 | heat(i,k)=0. |
|---|
| 114 | cool(i,k)=0. |
|---|
| 115 | ENDDO |
|---|
| 116 | ENDDO |
|---|
| 117 | |
|---|
| 118 | c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++ |
|---|
| 119 | DO j = 1, klon |
|---|
| 120 | |
|---|
| 121 | c====================================================================== |
|---|
| 122 | c Initialisations |
|---|
| 123 | c --------------- |
|---|
| 124 | |
|---|
| 125 | DO k = 1, klev |
|---|
| 126 | zheat(k) = 0.0 |
|---|
| 127 | zcool(k) = 0.0 |
|---|
| 128 | ENDDO |
|---|
| 129 | DO k = 1, klev+1 |
|---|
| 130 | ZFLNET(k) = 0.0 |
|---|
| 131 | ZFSNET(k) = 0.0 |
|---|
| 132 | ENDDO |
|---|
| 133 | ztopsw = 0.0 |
|---|
| 134 | ztoplw = 0.0 |
|---|
| 135 | zsolsw = 0.0 |
|---|
| 136 | zsollw = 0.0 |
|---|
| 137 | zsollwdown = 0.0 |
|---|
| 138 | |
|---|
| 139 | zfract = fract(j) |
|---|
| 140 | zrmu0 = rmu0(j) |
|---|
| 141 | |
|---|
| 142 | DO k = 1, klev+1 |
|---|
| 143 | PPB(k) = paprs(j,k)/1.e5 |
|---|
| 144 | ENDDO |
|---|
| 145 | |
|---|
| 146 | pt0(0) = tsol(j) |
|---|
| 147 | DO k = 1, klev |
|---|
| 148 | pt0(k) = t(j,k) |
|---|
| 149 | ENDDO |
|---|
| 150 | pt0(klev+1) = 0. |
|---|
| 151 | |
|---|
| 152 | DO k = 0,klev+1 |
|---|
| 153 | DO i = 0,klev+1 |
|---|
| 154 | psi(i,k) = 0. ! positif quand nrj de i->k |
|---|
| 155 | deltapsi(i,k) = 0. |
|---|
| 156 | ENDDO |
|---|
| 157 | ENDDO |
|---|
| 158 | |
|---|
| 159 | c====================================================================== |
|---|
| 160 | c Getting psi and deltapsi |
|---|
| 161 | c ------------------------ |
|---|
| 162 | |
|---|
| 163 | c Planck function |
|---|
| 164 | c --------------- |
|---|
| 165 | do band=1,nnuve |
|---|
| 166 | do k=0,klev |
|---|
| 167 | c B(T,l) = al/(exp(bl/T)-1) |
|---|
| 168 | y(k,band) = exp(bl(band)/pt0(k))-1. |
|---|
| 169 | bplck(k,band) = al(band)/(y(k,band)) |
|---|
| 170 | zdblay(k,band)= al(band)*bl(band)*exp(bl(band)/pt0(k))/ |
|---|
| 171 | . ((pt0(k)*pt0(k))*(y(k,band)*y(k,band))) |
|---|
| 172 | enddo |
|---|
| 173 | bplck(klev+1,band) = 0.0 |
|---|
| 174 | zdblay(klev+1,band)= 0.0 |
|---|
| 175 | enddo |
|---|
| 176 | |
|---|
| 177 | c finding the right matrixes |
|---|
| 178 | c -------------------------- |
|---|
| 179 | mat0 = 0 |
|---|
| 180 | if (nlatve.eq.1) then ! clouds are taken as uniform |
|---|
| 181 | lat=1 |
|---|
| 182 | else |
|---|
| 183 | if (abs(rlatd(j)).le.50.) then |
|---|
| 184 | lat=1 |
|---|
| 185 | elseif (abs(rlatd(j)).le.60.) then |
|---|
| 186 | lat=2 |
|---|
| 187 | elseif (abs(rlatd(j)).le.70.) then |
|---|
| 188 | lat=3 |
|---|
| 189 | elseif (abs(rlatd(j)).le.80.) then |
|---|
| 190 | lat=4 |
|---|
| 191 | else |
|---|
| 192 | lat=5 |
|---|
| 193 | endif |
|---|
| 194 | endif |
|---|
| 195 | |
|---|
| 196 | ips0=0 |
|---|
| 197 | do ips=1,nbpsve(lat)-1 |
|---|
| 198 | if ( (psurfve(ips,lat).ge.paprs(j,1)) |
|---|
| 199 | . .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then |
|---|
| 200 | ips0 = ips |
|---|
| 201 | c print*,'ig=',j,' ips0=',ips |
|---|
| 202 | factp = (paprs(j,1) -psurfve(ips0,lat)) |
|---|
| 203 | . /(psurfve(ips0+1,lat)-psurfve(ips0,lat)) |
|---|
| 204 | exit |
|---|
| 205 | endif |
|---|
| 206 | enddo |
|---|
| 207 | isza0=0 |
|---|
| 208 | if (nbszave(lat).gt.1) then |
|---|
| 209 | do isza=1,nbszave(lat)-1 |
|---|
| 210 | if ( (szave(isza,lat).ge.zrmu0) |
|---|
| 211 | . .and.(szave(isza+1,lat).lt.zrmu0) ) then |
|---|
| 212 | isza0 = isza |
|---|
| 213 | c print*,'ig=',j,' isza0=',isza |
|---|
| 214 | factz = (zrmu0 -szave(isza0,lat)) |
|---|
| 215 | . /(szave(isza0+1,lat)-szave(isza0,lat)) |
|---|
| 216 | exit |
|---|
| 217 | endif |
|---|
| 218 | enddo |
|---|
| 219 | else ! Only one sza, no interpolation |
|---|
| 220 | isza0=-99 |
|---|
| 221 | endif |
|---|
| 222 | |
|---|
| 223 | if ((ips0.eq.0).or.(isza0.eq.0)) then |
|---|
| 224 | write(*,*) 'Finding the right matrix in radlwsw' |
|---|
| 225 | print*,'Interpolation problem, grid point ig=',j |
|---|
| 226 | print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0 |
|---|
| 227 | stop |
|---|
| 228 | endif |
|---|
| 229 | |
|---|
| 230 | if (isza0.eq.-99) then |
|---|
| 231 | mat0 = indexve(lat)+ips0 |
|---|
| 232 | else |
|---|
| 233 | mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0 |
|---|
| 234 | endif |
|---|
| 235 | |
|---|
| 236 | c interpolation of ksi and computation of psi,deltapsi |
|---|
| 237 | c ---------------------------------------------------- |
|---|
| 238 | do band=1,nnuve |
|---|
| 239 | do k=0,klev+1 |
|---|
| 240 | do i=0,klev+1 |
|---|
| 241 | if (isza0.eq.-99) then |
|---|
| 242 | ksi = ksive(i,k,band,mat0)*(1-factp) |
|---|
| 243 | . +ksive(i,k,band,mat0+1)*factp |
|---|
| 244 | else |
|---|
| 245 | ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz) |
|---|
| 246 | . +ksive(i,k,band,mat0+1)*factp *(1-factz) |
|---|
| 247 | . +ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz |
|---|
| 248 | . +ksive(i,k,band,mat0+nbpsve(lat)+1)*factp *factz |
|---|
| 249 | endif |
|---|
| 250 | |
|---|
| 251 | psi(i,k) = psi(i,k) + |
|---|
| 252 | . RPI*ksi*(bplck(i,band)-bplck(k,band)) |
|---|
| 253 | deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band) |
|---|
| 254 | enddo |
|---|
| 255 | enddo |
|---|
| 256 | enddo |
|---|
| 257 | |
|---|
| 258 | c====================================================================== |
|---|
| 259 | c LW call |
|---|
| 260 | c--------- |
|---|
| 261 | temp(1:klev)=t(j,1:klev) |
|---|
| 262 | CALL LW_venus_ve( |
|---|
| 263 | . PPB,temp,psi,deltapsi, |
|---|
| 264 | . zcool, |
|---|
| 265 | . ztoplw,zsollw, |
|---|
| 266 | . zsollwdown,ZFLNET) |
|---|
| 267 | |
|---|
| 268 | c--------- |
|---|
| 269 | c SW call |
|---|
| 270 | c--------- |
|---|
| 271 | znivs=zzlev(j,:) |
|---|
| 272 | c CALL SW_venus_ve(zrmu0, zfract, |
|---|
| 273 | c S PPB,temp,znivs, |
|---|
| 274 | c S zheat, |
|---|
| 275 | c S ztopsw,zsolsw,ZFSNET) |
|---|
| 276 | |
|---|
| 277 | CALL SW_venus_dc(zrmu0, zfract, |
|---|
| 278 | S PPB,temp, |
|---|
| 279 | S zheat, |
|---|
| 280 | S ztopsw,zsolsw,ZFSNET) |
|---|
| 281 | |
|---|
| 282 | c====================================================================== |
|---|
| 283 | radsol(j) = zsolsw - zsollw ! + vers bas |
|---|
| 284 | topsw(j) = ztopsw ! + vers bas |
|---|
| 285 | toplw(j) = ztoplw ! + vers haut |
|---|
| 286 | solsw(j) = zsolsw ! + vers bas |
|---|
| 287 | sollw(j) = -zsollw ! + vers bas |
|---|
| 288 | sollwdown(j) = zsollwdown ! + vers bas |
|---|
| 289 | |
|---|
| 290 | DO k = 1, klev+1 |
|---|
| 291 | lwnet (j,k) = ZFLNET(k) |
|---|
| 292 | swnet (j,k) = ZFSNET(k) |
|---|
| 293 | ENDDO |
|---|
| 294 | |
|---|
| 295 | DO k = 1, klev |
|---|
| 296 | heat (j,k) = zheat(k) |
|---|
| 297 | cool (j,k) = zcool(k) |
|---|
| 298 | ENDDO |
|---|
| 299 | c |
|---|
| 300 | ENDDO ! of DO j = 1, klon |
|---|
| 301 | c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++ |
|---|
| 302 | ! for tests: write output fields... |
|---|
| 303 | ! call writefield_phy('radlwsw_heat',heat,klev) |
|---|
| 304 | ! call writefield_phy('radlwsw_cool',cool,klev) |
|---|
| 305 | ! call writefield_phy('radlwsw_radsol',radsol,1) |
|---|
| 306 | ! call writefield_phy('radlwsw_topsw',topsw,1) |
|---|
| 307 | ! call writefield_phy('radlwsw_toplw',toplw,1) |
|---|
| 308 | ! call writefield_phy('radlwsw_solsw',solsw,1) |
|---|
| 309 | ! call writefield_phy('radlwsw_sollw',sollw,1) |
|---|
| 310 | ! call writefield_phy('radlwsw_sollwdown',sollwdown,1) |
|---|
| 311 | ! call writefield_phy('radlwsw_swnet',swnet,klev+1) |
|---|
| 312 | ! call writefield_phy('radlwsw_lwnet',lwnet,klev+1) |
|---|
| 313 | |
|---|
| 314 | c tests |
|---|
| 315 | |
|---|
| 316 | c j = klon/2 |
|---|
| 317 | c j = 1 |
|---|
| 318 | c print*,'mu0=',rmu0(j) |
|---|
| 319 | c print*,' net flux vis HEAT(K/Eday)' |
|---|
| 320 | c do k=1,klev |
|---|
| 321 | c print*,k,ZFSNET(k),heat(j,k)*86400. |
|---|
| 322 | c enddo |
|---|
| 323 | c print*,' net flux IR COOL(K/Eday)' |
|---|
| 324 | c do k=1,klev |
|---|
| 325 | c print*,k,ZFLNET(k),cool(j,k)*86400. |
|---|
| 326 | c enddo |
|---|
| 327 | |
|---|
| 328 | firstcall = .false. |
|---|
| 329 | RETURN |
|---|
| 330 | END |
|---|
| 331 | |
|---|
