source: trunk/libf/phyvenus/ballon.F @ 91

Last change on this file since 91 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

  • libf/phytitan : physique de Titan
  • libf/chimtitan: chimie de Titan
  • libf/phyvenus : physique de Venus
File size: 12.3 KB
Line 
1      subroutine ballon (iam,dtphys,rjour,rsec,plat,plon,
2     i                   temp, p, u, v, geop)
3
4      implicit none
5
6c======================================================================
7c Auteur: S. Lebonnois (LMD/CNRS) date: 20091201
8c Object: Compute balloon trajectories.
9C No outputs, every quantities are written on the iam+ Files.
10c
11c Called by physiq.F if flag ballons activated:
12c
13c      integer ballons
14c (...)
15c      ballons  = 1         ! (in initialisations)
16c (...)
17C  OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
18C  DES BALLONS
19c      if (ballons.eq.1) then
20c      open(30,file='ballons-lat.out',form='formatted')
21c      open(31,file='ballons-lon.out',form='formatted')
22c      open(32,file='ballons-u.out',form='formatted')
23c      open(33,file='ballons-v.out',form='formatted')
24c      open(34,file='ballons-alt.out',form='formatted')
25c      write(*,*)'Ouverture des ballons*.out'
26c      endif !ballons
27c (...)
28C        CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
29CC    C               t,pplay,u,v,pphi)   ! alt above surface (smoothed for GCM)
30C     C               t,pplay,u,v,zphi)   ! alt above planet average radius
31c (...)
32C  FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
33C  DES BALLONS
34c      if (ballons.eq.1) then
35c        write(*,*)'Fermeture des ballons*.out'
36c        close(30)                                     
37c        close(31)                                     
38c        close(32)                                     
39c        close(33)                                     
40c        close(34)                                     
41c      endif !ballons
42c
43C======================================================================
44c Explicit Arguments:
45c ==================
46c iam-----input-I-File number where latitudes are written
47c                 It is a formatted file that has been opened
48c                 in physiq.F
49c   other files: iam+1=longitudes
50c                iam+2=zonal speeds
51c                iam+3=meridional speeds
52c                iam+4=altitudes
53c dtphys--input-R-pas de temps physique
54c rjour---input-R-Jour compte depuis le debut de la simu (run.def)
55c rsec----input-R-Seconde de la journee
56c plat ---input-R-Latitude en degres
57c plon ---input-R-Longitude en degres
58c temp----input-R-Temperature (K) at model levels
59c p-------input-R-Pressure (Pa) at model levels
60c u-------input-R-Horizontal wind (m/s)
61c v-------input-R-Meridional wind (m/s)
62c geop----input-R-Geopotential !! above surface OR average radius
63c
64c
65c Implicit Arguments:
66c ===================
67c
68c iim--common-I: Number of longitude intervals
69c jjm--common-I: Number of latitude intervals
70c klon-common-I: Number of points seen by the physics
71c                iim*(jjm-1)+2 for instance
72c klev-common-I: Number of vertical layers
73c RPI,RKBOL--common-R: Pi, KBoltzman
74c RDAY,RA,RG-common-R: day length in s, planet radius, gravity
75c======================================================================
76c Local Variables:
77c ================
78c
79c nb    ---I: number of balloons (parameter)
80c phib  ---R: Latitude of balloon in radians
81c lamb  ---R: Longitude of balloon in radians
82c lognb ---R: log(density) of balloon
83c ub    ---R: zonal speed of balloon
84c vb    ---R: meridional speed of balloon
85c altb  ---R: altitude of balloon
86c zlat  ---R: Latitude in radians
87c zlon  ---R: Longitude in radians
88c logn  ---R: log(density)
89c alt   ---R: altitude !! above surface OR average radius
90c ull   ---R: zonal wind for one balloon on the lognb surface
91c vll   ---R: meridional wind for one balloon on the lognb surface
92c aal   ---R: altitude for one balloon on the lognb surface
93c======================================================================
94
95#include "dimensions.h"
96#include "dimphy.h"
97#include "YOMCST.h"
98c
99c ARGUMENTS
100c
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)
105c
106c Variables locales:
107c
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
123c 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
136c---------------------------------------------
137C INITIALIZATIONS
138c---------------------------------------------
139      if (first) then
140
141      print*,"BALLOONS ACTIVATED"
142
143C 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
151C 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
158c verif init     lat de 90 à -90, lon de -180 à 180
159c     print*,"Latitudes:",zlat*180./RPI
160c     print*,"Longitudes:",zlon*180./RPI
161c     stop
162
163c 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
174c     print*,"log(density) profile:"
175c     do l=1,klev
176c        print*,logn(klon/2,l)
177c     enddo
178c     stop !verif init
179
180      first=.false.
181      endif ! first
182c---------------------------------------------
183
184c-------------------------------------------------
185c loop over the balloons
186c-------------------------------------------------
187      do n=1,nb
188
189c Interpolation in altitudes
190c-------------------------------------------------
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
202c Interpolation in latitudes and longitudes
203c-------------------------------------------
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
208c-------------------------------------------------
209
210c-------------------------------------------------
211c Output of positions and speed at time
212c-------------------------------------------------
213
214c 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
218c Venus regardee a l'envers: v inversee
219      write(iam+3,100) time, (-1)*vb
220      write(iam+4,100) time, altb
221c     stop !verif init
222
223c !!!!!!!!!!!!!!!! nb !!!!!!!!!!!!!!!!!
224100   format(E14.7,20(1x,E12.5))
225
226c-------------------------------------------------
227c Implementation: positions at time+dt
228c RK order Nrk
229c-------------------------------------------------
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
251c======================================================================
252c======================================================================
253c======================================================================
254
255      subroutine wind_interp(map_u,map_v,map_a,latit,longit,
256     .                       phi,lam,ubal,vbal,abal)
257
258      implicit none
259
260c======================================================================
261c Auteur: S. Lebonnois (LMD/CNRS) date: 20091201
262c Object: interpolate balloon speed from its position.
263C======================================================================
264c Explicit Arguments:
265c ==================
266c map_u ---R: zonal wind on the lognb surface
267c map_v ---R: meridional wind on the lognb surface
268c map_a ---R: altitude on the lognb surface
269c latit ---R: Latitude in radians
270c longit---R: Longitude in radians
271c phi   ---R: Latitude of balloon in radians
272c lam   ---R: Longitude of balloon in radians
273c ubal  ---R: zonal speed of balloon
274c vbal  ---R: meridional speed of balloon
275c abal  ---R: altitude of balloon
276c======================================================================
277c Local Variables:
278c ================
279c
280c ujj   ---R: zonal wind interpolated in latitude
281c vjj   ---R: meridional wind interpolated in latitude
282c ajj   ---R: altitude interpolated in latitude
283c======================================================================
284
285#include "dimensions.h"
286#include "dimphy.h"
287#include "YOMCST.h"
288c
289c ARGUMENTS
290c
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
294c
295c Variables locales:
296c
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
302c Interpolation in latitudes
303c-------------------------------------------------
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
310c 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
317c 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
325c 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
338c Interpolation in longitudes
339c-------------------------------------------------
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
351c======================================================================
352c======================================================================
353c======================================================================
354
355      subroutine pos_implem(phi,lam,ubal,vbal,dt)
356
357      implicit none
358
359c======================================================================
360c Auteur: S. Lebonnois (LMD/CNRS) date: 20091201
361c Object: implementation of balloon position.
362C======================================================================
363c Explicit Arguments:
364c ==================
365c phi   ---R: Latitude of balloon in radians
366c lam   ---R: Longitude of balloon in radians
367c ubal  ---R: zonal speed of balloon
368c vbal  ---R: meridional speed of balloon
369c dt    ---R: time step
370c======================================================================
371
372#include "dimensions.h"
373#include "dimphy.h"
374#include "YOMCST.h"
375c
376c ARGUMENTS
377c
378      real phi,lam,ubal,vbal,abal,dt
379
380c incrementation longitude
381        lam = lam + ubal*dt/(RA*cos(phi))
382c 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
389c incrementation latitude
390        phi = phi + vbal*dt/RA
391c 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
Note: See TracBrowser for help on using the repository browser.