| 1 | subroutine ballon (iam,dtphys,rjour,rsec,plat,plon, |
|---|
| 2 | i temp, p, u, v, geop) |
|---|
| 3 | |
|---|
| 4 | use dimphy |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | c====================================================================== |
|---|
| 8 | c Auteur: S. Lebonnois (LMD/CNRS) date: 20091201 |
|---|
| 9 | c Object: Compute balloon trajectories. |
|---|
| 10 | C No outputs, every quantities are written on the iam+ Files. |
|---|
| 11 | c |
|---|
| 12 | c Called by physiq.F if flag ballons activated: |
|---|
| 13 | c |
|---|
| 14 | c integer ballons |
|---|
| 15 | c (...) |
|---|
| 16 | c ballons = 1 ! (in initialisations) |
|---|
| 17 | c (...) |
|---|
| 18 | C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES |
|---|
| 19 | C DES BALLONS |
|---|
| 20 | c if (ballons.eq.1) then |
|---|
| 21 | c open(30,file='ballons-lat.out',form='formatted') |
|---|
| 22 | c open(31,file='ballons-lon.out',form='formatted') |
|---|
| 23 | c open(32,file='ballons-u.out',form='formatted') |
|---|
| 24 | c open(33,file='ballons-v.out',form='formatted') |
|---|
| 25 | c open(34,file='ballons-alt.out',form='formatted') |
|---|
| 26 | c write(*,*)'Ouverture des ballons*.out' |
|---|
| 27 | c endif !ballons |
|---|
| 28 | c (...) |
|---|
| 29 | C CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond, |
|---|
| 30 | CC C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM) |
|---|
| 31 | C C t,pplay,u,v,zphi) ! alt above planet average radius |
|---|
| 32 | c (...) |
|---|
| 33 | C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES |
|---|
| 34 | C DES BALLONS |
|---|
| 35 | c if (ballons.eq.1) then |
|---|
| 36 | c write(*,*)'Fermeture des ballons*.out' |
|---|
| 37 | c close(30) |
|---|
| 38 | c close(31) |
|---|
| 39 | c close(32) |
|---|
| 40 | c close(33) |
|---|
| 41 | c close(34) |
|---|
| 42 | c endif !ballons |
|---|
| 43 | c |
|---|
| 44 | C====================================================================== |
|---|
| 45 | c Explicit Arguments: |
|---|
| 46 | c ================== |
|---|
| 47 | c iam-----input-I-File number where latitudes are written |
|---|
| 48 | c It is a formatted file that has been opened |
|---|
| 49 | c in physiq.F |
|---|
| 50 | c other files: iam+1=longitudes |
|---|
| 51 | c iam+2=zonal speeds |
|---|
| 52 | c iam+3=meridional speeds |
|---|
| 53 | c iam+4=altitudes |
|---|
| 54 | c dtphys--input-R-pas de temps physique |
|---|
| 55 | c rjour---input-R-Jour compte depuis le debut de la simu (run.def) |
|---|
| 56 | c rsec----input-R-Seconde de la journee |
|---|
| 57 | c plat ---input-R-Latitude en degres |
|---|
| 58 | c plon ---input-R-Longitude en degres |
|---|
| 59 | c temp----input-R-Temperature (K) at model levels |
|---|
| 60 | c p-------input-R-Pressure (Pa) at model levels |
|---|
| 61 | c u-------input-R-Horizontal wind (m/s) |
|---|
| 62 | c v-------input-R-Meridional wind (m/s) |
|---|
| 63 | c geop----input-R-Geopotential !! above surface OR average radius |
|---|
| 64 | c |
|---|
| 65 | c |
|---|
| 66 | c Implicit Arguments: |
|---|
| 67 | c =================== |
|---|
| 68 | c |
|---|
| 69 | c iim--common-I: Number of longitude intervals |
|---|
| 70 | c jjm--common-I: Number of latitude intervals |
|---|
| 71 | c klon-common-I: Number of points seen by the physics |
|---|
| 72 | c iim*(jjm-1)+2 for instance |
|---|
| 73 | c klev-common-I: Number of vertical layers |
|---|
| 74 | c RPI,RKBOL--common-R: Pi, KBoltzman |
|---|
| 75 | c RDAY,RA,RG-common-R: day length in s, planet radius, gravity |
|---|
| 76 | c====================================================================== |
|---|
| 77 | c Local Variables: |
|---|
| 78 | c ================ |
|---|
| 79 | c |
|---|
| 80 | c nb ---I: number of balloons (parameter) |
|---|
| 81 | c phib ---R: Latitude of balloon in radians |
|---|
| 82 | c lamb ---R: Longitude of balloon in radians |
|---|
| 83 | c lognb ---R: log(density) of balloon |
|---|
| 84 | c ub ---R: zonal speed of balloon |
|---|
| 85 | c vb ---R: meridional speed of balloon |
|---|
| 86 | c altb ---R: altitude of balloon |
|---|
| 87 | c zlat ---R: Latitude in radians |
|---|
| 88 | c zlon ---R: Longitude in radians |
|---|
| 89 | c logn ---R: log(density) |
|---|
| 90 | c alt ---R: altitude !! above surface OR average radius |
|---|
| 91 | c ull ---R: zonal wind for one balloon on the lognb surface |
|---|
| 92 | c vll ---R: meridional wind for one balloon on the lognb surface |
|---|
| 93 | c aal ---R: altitude for one balloon on the lognb surface |
|---|
| 94 | c====================================================================== |
|---|
| 95 | |
|---|
| 96 | #include "dimensions.h" |
|---|
| 97 | #include "YOMCST.h" |
|---|
| 98 | c |
|---|
| 99 | c ARGUMENTS |
|---|
| 100 | c |
|---|
| 101 | INTEGER iam |
|---|
| 102 | REAL dtphys,rjour,rsec,plat(klon),plon(klon) |
|---|
| 103 | REAL temp(klon,klev),p(klon,klev) |
|---|
| 104 | REAL u(klon,klev),v(klon,klev),geop(klon,klev) |
|---|
| 105 | c |
|---|
| 106 | c Variables locales: |
|---|
| 107 | c |
|---|
| 108 | INTEGER i,j,k,l,nb,n |
|---|
| 109 | parameter (nb=20) !! Adjust the format on line 100 !! |
|---|
| 110 | INTEGER jj,ii,ll |
|---|
| 111 | |
|---|
| 112 | REAL zlon(iim+1),zlat(jjm+1) |
|---|
| 113 | save zlon,zlat |
|---|
| 114 | |
|---|
| 115 | REAL time |
|---|
| 116 | REAL logn(klon,klev),ull(klon),vll(klon) |
|---|
| 117 | REAL alt(klon,klev),aal(klon) |
|---|
| 118 | real ub(nb),vb(nb),phib(nb),lamb(nb),lognb(nb),altb(nb) |
|---|
| 119 | save phib,lamb,lognb |
|---|
| 120 | |
|---|
| 121 | REAL factalt |
|---|
| 122 | |
|---|
| 123 | c RungeKutta order - If not RK, Nrk=1 |
|---|
| 124 | integer Nrk,irk |
|---|
| 125 | parameter (Nrk=1) |
|---|
| 126 | real dtrk |
|---|
| 127 | |
|---|
| 128 | logical first |
|---|
| 129 | save first |
|---|
| 130 | data first/.true./ |
|---|
| 131 | |
|---|
| 132 | time = rjour*RDAY+rsec |
|---|
| 133 | logn(:,:) = log10(p(:,:)/(RKBOL*temp(:,:))) |
|---|
| 134 | alt(:,:) = geop(:,:)/RG |
|---|
| 135 | |
|---|
| 136 | c--------------------------------------------- |
|---|
| 137 | C INITIALIZATIONS |
|---|
| 138 | c--------------------------------------------- |
|---|
| 139 | if (first) then |
|---|
| 140 | |
|---|
| 141 | print*,"BALLOONS ACTIVATED" |
|---|
| 142 | |
|---|
| 143 | C Latitudes: |
|---|
| 144 | zlat(1)=plat(1)*RPI/180. |
|---|
| 145 | do j = 2,jjm |
|---|
| 146 | k=(j-2)*iim+2 |
|---|
| 147 | zlat(j)=plat(k)*RPI/180. |
|---|
| 148 | enddo |
|---|
| 149 | zlat(jjm+1)=plat(klon)*RPI/180. |
|---|
| 150 | |
|---|
| 151 | C Longitudes: |
|---|
| 152 | do i = 1,iim |
|---|
| 153 | k=i+1 |
|---|
| 154 | zlon(i)=plon(k)*RPI/180. |
|---|
| 155 | enddo |
|---|
| 156 | zlon(iim+1)=zlon(1)+2.*RPI |
|---|
| 157 | |
|---|
| 158 | c verif init lat de 90 à -90, lon de -180 à 180 |
|---|
| 159 | c print*,"Latitudes:",zlat*180./RPI |
|---|
| 160 | c print*,"Longitudes:",zlon*180./RPI |
|---|
| 161 | c stop |
|---|
| 162 | |
|---|
| 163 | c initial positions of balloons (in degrees for lat/lon) |
|---|
| 164 | do j=1,5 |
|---|
| 165 | do i=1,4 |
|---|
| 166 | k=(j-1)*4+i |
|---|
| 167 | phib(k)= (j-1)*20.*RPI/180. |
|---|
| 168 | lamb(k)= (i-3)*90.*RPI/180. ! de -180 à 90 |
|---|
| 169 | lognb(k)= log10(5.e4/(RKBOL*300.)) ! ~55km in VIRA model |
|---|
| 170 | enddo |
|---|
| 171 | enddo |
|---|
| 172 | print*,"Balloon density (m^-3)=",10.**(lognb(1)) |
|---|
| 173 | |
|---|
| 174 | c print*,"log(density) profile:" |
|---|
| 175 | c do l=1,klev |
|---|
| 176 | c print*,logn(klon/2,l) |
|---|
| 177 | c enddo |
|---|
| 178 | c stop !verif init |
|---|
| 179 | |
|---|
| 180 | first=.false. |
|---|
| 181 | endif ! first |
|---|
| 182 | c--------------------------------------------- |
|---|
| 183 | |
|---|
| 184 | c------------------------------------------------- |
|---|
| 185 | c loop over the balloons |
|---|
| 186 | c------------------------------------------------- |
|---|
| 187 | do n=1,nb |
|---|
| 188 | |
|---|
| 189 | c Interpolation in altitudes |
|---|
| 190 | c------------------------------------------------- |
|---|
| 191 | do k=1,klon |
|---|
| 192 | ll=1 ! en bas |
|---|
| 193 | do l=2,klev |
|---|
| 194 | if (lognb(n).lt.logn(k,l)) ll=l |
|---|
| 195 | enddo |
|---|
| 196 | factalt= (lognb(n)-logn(k,ll))/(logn(k,ll+1)-logn(k,ll)) |
|---|
| 197 | ull(k) = u(k,ll+1)*factalt + u(k,ll)*(1-factalt) |
|---|
| 198 | vll(k) = v(k,ll+1)*factalt + v(k,ll)*(1-factalt) |
|---|
| 199 | aal(k) = alt(k,ll+1)*factalt + alt(k,ll)*(1-factalt) |
|---|
| 200 | enddo |
|---|
| 201 | |
|---|
| 202 | c Interpolation in latitudes and longitudes |
|---|
| 203 | c------------------------------------------- |
|---|
| 204 | call wind_interp(ull,vll,aal,zlat,zlon, |
|---|
| 205 | . phib(n),lamb(n),ub(n),vb(n),altb(n)) |
|---|
| 206 | |
|---|
| 207 | enddo ! over balloons |
|---|
| 208 | c------------------------------------------------- |
|---|
| 209 | |
|---|
| 210 | c------------------------------------------------- |
|---|
| 211 | c Output of positions and speed at time |
|---|
| 212 | c------------------------------------------------- |
|---|
| 213 | |
|---|
| 214 | c Venus regardee a l'envers: lon et lat inversees |
|---|
| 215 | write(iam, 100) time, (-1)*phib*180./RPI |
|---|
| 216 | write(iam+1,100) time, (-1)*lamb*180./RPI |
|---|
| 217 | write(iam+2,100) time, ub |
|---|
| 218 | c Venus regardee a l'envers: v inversee |
|---|
| 219 | write(iam+3,100) time, (-1)*vb |
|---|
| 220 | write(iam+4,100) time, altb |
|---|
| 221 | c stop !verif init |
|---|
| 222 | |
|---|
| 223 | c !!!!!!!!!!!!!!!! nb !!!!!!!!!!!!!!!!! |
|---|
| 224 | 100 format(E14.7,20(1x,E12.5)) |
|---|
| 225 | |
|---|
| 226 | c------------------------------------------------- |
|---|
| 227 | c Implementation: positions at time+dt |
|---|
| 228 | c RK order Nrk |
|---|
| 229 | c------------------------------------------------- |
|---|
| 230 | |
|---|
| 231 | dtrk = dtphys/Nrk |
|---|
| 232 | time=time+dtrk |
|---|
| 233 | |
|---|
| 234 | do n=1,nb |
|---|
| 235 | call pos_implem(phib(n),lamb(n),ub(n),vb(n),dtrk) |
|---|
| 236 | enddo |
|---|
| 237 | |
|---|
| 238 | if (Nrk.gt.1) then |
|---|
| 239 | do irk=2,Nrk |
|---|
| 240 | do n=1,nb |
|---|
| 241 | time=time+dtrk |
|---|
| 242 | call wind_interp(ull,vll,aal,zlat,zlon, |
|---|
| 243 | . phib(n),lamb(n),ub(n),vb(n),altb(n)) |
|---|
| 244 | call pos_implem(phib(n),lamb(n),ub(n),vb(n),dtrk) |
|---|
| 245 | enddo |
|---|
| 246 | enddo |
|---|
| 247 | endif |
|---|
| 248 | |
|---|
| 249 | end |
|---|
| 250 | |
|---|
| 251 | c====================================================================== |
|---|
| 252 | c====================================================================== |
|---|
| 253 | c====================================================================== |
|---|
| 254 | |
|---|
| 255 | subroutine wind_interp(map_u,map_v,map_a,latit,longit, |
|---|
| 256 | . phi,lam,ubal,vbal,abal) |
|---|
| 257 | |
|---|
| 258 | use dimphy |
|---|
| 259 | implicit none |
|---|
| 260 | |
|---|
| 261 | c====================================================================== |
|---|
| 262 | c Auteur: S. Lebonnois (LMD/CNRS) date: 20091201 |
|---|
| 263 | c Object: interpolate balloon speed from its position. |
|---|
| 264 | C====================================================================== |
|---|
| 265 | c Explicit Arguments: |
|---|
| 266 | c ================== |
|---|
| 267 | c map_u ---R: zonal wind on the lognb surface |
|---|
| 268 | c map_v ---R: meridional wind on the lognb surface |
|---|
| 269 | c map_a ---R: altitude on the lognb surface |
|---|
| 270 | c latit ---R: Latitude in radians |
|---|
| 271 | c longit---R: Longitude in radians |
|---|
| 272 | c phi ---R: Latitude of balloon in radians |
|---|
| 273 | c lam ---R: Longitude of balloon in radians |
|---|
| 274 | c ubal ---R: zonal speed of balloon |
|---|
| 275 | c vbal ---R: meridional speed of balloon |
|---|
| 276 | c abal ---R: altitude of balloon |
|---|
| 277 | c====================================================================== |
|---|
| 278 | c Local Variables: |
|---|
| 279 | c ================ |
|---|
| 280 | c |
|---|
| 281 | c ujj ---R: zonal wind interpolated in latitude |
|---|
| 282 | c vjj ---R: meridional wind interpolated in latitude |
|---|
| 283 | c ajj ---R: altitude interpolated in latitude |
|---|
| 284 | c====================================================================== |
|---|
| 285 | |
|---|
| 286 | #include "dimensions.h" |
|---|
| 287 | #include "YOMCST.h" |
|---|
| 288 | c |
|---|
| 289 | c ARGUMENTS |
|---|
| 290 | c |
|---|
| 291 | real map_u(klon),map_v(klon),map_a(klon) |
|---|
| 292 | real latit(jjm+1),longit(iim) |
|---|
| 293 | real phi,lam,ubal,vbal,abal |
|---|
| 294 | c |
|---|
| 295 | c Variables locales: |
|---|
| 296 | c |
|---|
| 297 | INTEGER i,j,k |
|---|
| 298 | INTEGER jj,ii |
|---|
| 299 | REAL ujj(iim+1),vjj(iim+1),ajj(iim+1) |
|---|
| 300 | REAL factlat,factlon |
|---|
| 301 | |
|---|
| 302 | c Interpolation in latitudes |
|---|
| 303 | c------------------------------------------------- |
|---|
| 304 | jj=1 ! POLE NORD |
|---|
| 305 | do j=2,jjm |
|---|
| 306 | if (phi.lt.latit(j)) jj=j |
|---|
| 307 | enddo |
|---|
| 308 | factlat = (phi-latit(jj))/(latit(jj+1)-latit(jj)) |
|---|
| 309 | |
|---|
| 310 | c pole nord |
|---|
| 311 | if (jj.eq.1) then |
|---|
| 312 | do i=1,iim |
|---|
| 313 | ujj(i) = map_u(i+1)*factlat + map_u(1)*(1-factlat) |
|---|
| 314 | vjj(i) = map_v(i+1)*factlat + map_v(1)*(1-factlat) |
|---|
| 315 | ajj(i) = map_a(i+1)*factlat + map_a(1)*(1-factlat) |
|---|
| 316 | enddo |
|---|
| 317 | c pole sud |
|---|
| 318 | elseif (jj.eq.jjm) then |
|---|
| 319 | do i=1,iim |
|---|
| 320 | k = (jj-2)*iim+1+i |
|---|
| 321 | ujj(i) = map_u(klon)*factlat + map_u(k)*(1-factlat) |
|---|
| 322 | vjj(i) = map_v(klon)*factlat + map_v(k)*(1-factlat) |
|---|
| 323 | ajj(i) = map_a(klon)*factlat + map_a(k)*(1-factlat) |
|---|
| 324 | enddo |
|---|
| 325 | c autres latitudes |
|---|
| 326 | else |
|---|
| 327 | do i=1,iim |
|---|
| 328 | k = (jj-2)*iim+1+i |
|---|
| 329 | ujj(i) = map_u(k+iim)*factlat + map_u(k)*(1-factlat) |
|---|
| 330 | vjj(i) = map_v(k+iim)*factlat + map_v(k)*(1-factlat) |
|---|
| 331 | ajj(i) = map_a(k+iim)*factlat + map_a(k)*(1-factlat) |
|---|
| 332 | enddo |
|---|
| 333 | endif |
|---|
| 334 | ujj(iim+1)=ujj(1) |
|---|
| 335 | vjj(iim+1)=vjj(1) |
|---|
| 336 | ajj(iim+1)=ajj(1) |
|---|
| 337 | |
|---|
| 338 | c Interpolation in longitudes |
|---|
| 339 | c------------------------------------------------- |
|---|
| 340 | ii=1 ! lon=-180 |
|---|
| 341 | do i=2,iim |
|---|
| 342 | if (lam.gt.longit(i)) ii=i |
|---|
| 343 | enddo |
|---|
| 344 | factlon = (lam-longit(ii))/(longit(ii+1)-longit(ii)) |
|---|
| 345 | ubal = ujj(ii+1)*factlon + ujj(ii)*(1-factlon) |
|---|
| 346 | vbal = vjj(ii+1)*factlon + vjj(ii)*(1-factlon) |
|---|
| 347 | abal = ajj(ii+1)*factlon + ajj(ii)*(1-factlon) |
|---|
| 348 | |
|---|
| 349 | end |
|---|
| 350 | |
|---|
| 351 | c====================================================================== |
|---|
| 352 | c====================================================================== |
|---|
| 353 | c====================================================================== |
|---|
| 354 | |
|---|
| 355 | subroutine pos_implem(phi,lam,ubal,vbal,dt) |
|---|
| 356 | |
|---|
| 357 | use dimphy |
|---|
| 358 | implicit none |
|---|
| 359 | |
|---|
| 360 | c====================================================================== |
|---|
| 361 | c Auteur: S. Lebonnois (LMD/CNRS) date: 20091201 |
|---|
| 362 | c Object: implementation of balloon position. |
|---|
| 363 | C====================================================================== |
|---|
| 364 | c Explicit Arguments: |
|---|
| 365 | c ================== |
|---|
| 366 | c phi ---R: Latitude of balloon in radians |
|---|
| 367 | c lam ---R: Longitude of balloon in radians |
|---|
| 368 | c ubal ---R: zonal speed of balloon |
|---|
| 369 | c vbal ---R: meridional speed of balloon |
|---|
| 370 | c dt ---R: time step |
|---|
| 371 | c====================================================================== |
|---|
| 372 | |
|---|
| 373 | #include "dimensions.h" |
|---|
| 374 | #include "YOMCST.h" |
|---|
| 375 | c |
|---|
| 376 | c ARGUMENTS |
|---|
| 377 | c |
|---|
| 378 | real phi,lam,ubal,vbal,abal,dt |
|---|
| 379 | |
|---|
| 380 | c incrementation longitude |
|---|
| 381 | lam = lam + ubal*dt/(RA*cos(phi)) |
|---|
| 382 | c maintenue entre -PI et PI: |
|---|
| 383 | do while (lam.ge.RPI) |
|---|
| 384 | lam=lam-2*RPI |
|---|
| 385 | enddo |
|---|
| 386 | do while (lam.lt.(-1.*RPI)) |
|---|
| 387 | lam=lam+2*RPI |
|---|
| 388 | enddo |
|---|
| 389 | c incrementation latitude |
|---|
| 390 | phi = phi + vbal*dt/RA |
|---|
| 391 | c maintenue entre -PI/2 et PI/2: |
|---|
| 392 | if (phi.ge.( 0.5*RPI)) phi= RPI-phi |
|---|
| 393 | if (phi.le.(-0.5*RPI)) phi=-1.*RPI-phi |
|---|
| 394 | |
|---|
| 395 | end |
|---|
