source: trunk/LMDZ.VENUS/libf/phyvenus/ballon.F @ 2187

Last change on this file since 2187 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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