My Project
 All Classes Files Functions Variables Macros
advtrac.F90
Go to the documentation of this file.
1 ! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z lguez $
2 
3 SUBROUTINE advtrac(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk)
4  ! Auteur : F. Hourdin
5  !
6  ! Modif. P. Le Van (20/12/97)
7  ! F. Codron (10/99)
8  ! D. Le Croller (07/2001)
9  ! M.A Filiberti (04/2002)
10  !
11  USE infotrac
12  USE control_mod
13 
14 
15  IMPLICIT NONE
16  !
17  include "dimensions.h"
18  include "paramet.h"
19  include "comconst.h"
20  include "comvert.h"
21  include "comdissip.h"
22  include "comgeom2.h"
23  include "logic.h"
24  include "temps.h"
25  include "ener.h"
26  include "description.h"
27  include "iniprint.h"
28 
29  !-------------------------------------------------------------------
30  ! Arguments
31  !-------------------------------------------------------------------
32  ! Ajout PPM
33  !--------------------------------------------------------
34  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
35  !--------------------------------------------------------
36  INTEGER iapptrac
37  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
38  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
39  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
40  REAL pk(ip1jmp1,llm)
41  REAL flxw(ip1jmp1,llm)
42 
43  !-------------------------------------------------------------
44  ! Variables locales
45  !-------------------------------------------------------------
46 
47  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
48  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
49  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
50  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
51  INTEGER iadvtr
52  INTEGER ij,l,iq,iiq
53  REAL zdpmin, zdpmax
54  EXTERNAL minmax
55  SAVE iadvtr, massem, pbaruc, pbarvc
56  DATA iadvtr/0/
57  !----------------------------------------------------------
58  ! Rajouts pour PPM
59  !----------------------------------------------------------
60  INTEGER indice,n
61  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
62  REAL cflmaxz,aaa,bbb ! CFL maximum
63  REAL psppm(iim,jjp1) ! pression au sol
64  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
65  REAL qppm(iim*jjp1,llm,nqtot)
66  REAL fluxwppm(iim,jjp1,llm)
67  REAL apppm(llmp1), bpppm(llmp1)
68  LOGICAL dum,fill
69  DATA fill/.true./
70  DATA dum/.true./
71 
72  integer,save :: countcfl=0
73  real cflx(ip1jmp1,llm)
74  real cfly(ip1jm,llm)
75  real cflz(ip1jmp1,llm)
76  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
77 
78  IF(iadvtr.EQ.0) THEN
79  CALL initial0(ijp1llm,pbaruc)
80  CALL initial0(ijmllm,pbarvc)
81  ENDIF
82 
83  ! accumulation des flux de masse horizontaux
84  DO l=1,llm
85  DO ij = 1,ip1jmp1
86  pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
87  ENDDO
88  DO ij = 1,ip1jm
89  pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
90  ENDDO
91  ENDDO
92 
93  ! selection de la masse instantannee des mailles avant le transport.
94  IF(iadvtr.EQ.0) THEN
95 
96  CALL scopy(ip1jmp1*llm,masse,1,massem,1)
97  !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
98  !
99  ENDIF
100 
101  iadvtr = iadvtr+1
102  iapptrac = iadvtr
103 
104 
105  ! Test pour savoir si on advecte a ce pas de temps
106  IF ( iadvtr.EQ.iapp_tracvl ) THEN
107 
108  !c .. Modif P.Le Van ( 20/12/97 ) ....
109  !c
110 
111  ! traitement des flux de masse avant advection.
112  ! 1. calcul de w
113  ! 2. groupement des mailles pres du pole.
114 
115  CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
116 
117  ! ... Flux de masse diaganostiques traceurs
118  flxw = wg / REAL(iapp_tracvl)
119 
120  ! test sur l'eventuelle creation de valeurs negatives de la masse
121  DO l=1,llm-1
122  DO ij = iip2+1,ip1jm
123  zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
124  - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
125  + wg(ij,l+1) - wg(ij,l)
126  ENDDO
127  CALL scopy( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
128  DO ij = iip2,ip1jm
129  zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
130  ENDDO
131 
132 
133  CALL minmax( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
134 
135  IF(max(abs(zdpmin),abs(zdpmax)).GT.0.5) THEN
136  print*,'WARNING DP/P l=',l,' MIN:',zdpmin, &
137  ' MAX:', zdpmax
138  ENDIF
139 
140  ENDDO
141 
142 
143  !-------------------------------------------------------------------
144  ! Calcul des criteres CFL en X, Y et Z
145  !-------------------------------------------------------------------
146 
147  if (countcfl == 0. ) then
148  cflxmax(:)=0.
149  cflymax(:)=0.
150  cflzmax(:)=0.
151  endif
152 
153  countcfl=countcfl+iapp_tracvl
154  cflx(:,:)=0.
155  cfly(:,:)=0.
156  cflz(:,:)=0.
157  do l=1,llm
158  do ij=iip2,ip1jm-1
159  if (pbarug(ij,l)>=0.) then
160  cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
161  else
162  cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
163  endif
164  enddo
165  enddo
166  do l=1,llm
167  do ij=iip2,ip1jm-1,iip1
168  cflx(ij+iip1,l)=cflx(ij,l)
169  enddo
170  enddo
171 
172  do l=1,llm
173  do ij=1,ip1jm
174  if (pbarvg(ij,l)>=0.) then
175  cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
176  else
177  cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
178  endif
179  enddo
180  enddo
181 
182  do l=2,llm
183  do ij=1,ip1jm
184  if (wg(ij,l)>=0.) then
185  cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
186  else
187  cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
188  endif
189  enddo
190  enddo
191 
192  do l=1,llm
193  cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
194  cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
195  cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
196  enddo
197 
198 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
199  ! Par defaut, on sort le diagnostic des CFL tous les jours.
200  ! Si on veut le sortir a chaque pas d'advection en cas de plantage
201  ! if (countcfl==iapp_tracvl) then
202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203  if (countcfl==day_step) then
204  do l=1,llm
205  write(lunout,*) 'L, CFLmax ' &
206  ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
207  enddo
208  countcfl=0
209  endif
210 
211  !-------------------------------------------------------------------
212  ! Advection proprement dite (Modification Le Croller (07/2001)
213  !-------------------------------------------------------------------
214 
215  !----------------------------------------------------
216  ! Calcul des moyennes basées sur la masse
217  !----------------------------------------------------
218  call massbar(massem,massebx,masseby)
219 
220  !-----------------------------------------------------------
221  ! Appel des sous programmes d'advection
222  !-----------------------------------------------------------
223  do iq=1,nqtot
224  ! call clock(t_initial)
225  if(iadv(iq) == 0) cycle
226  ! ----------------------------------------------------------------
227  ! Schema de Van Leer I MUSCL
228  ! ----------------------------------------------------------------
229  if(iadv(iq).eq.10) THEN
230  call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
231 
232  ! ----------------------------------------------------------------
233  ! Schema "pseudo amont" + test sur humidite specifique
234  ! pour la vapeur d'eau. F. Codron
235  ! ----------------------------------------------------------------
236  else if(iadv(iq).eq.14) then
237  !
238  CALL vlspltqs( q(1,1,1), 2., massem, wg , &
239  pbarug,pbarvg,dtvr,p,pk,teta )
240  ! ----------------------------------------------------------------
241  ! Schema de Frederic Hourdin
242  ! ----------------------------------------------------------------
243  else if(iadv(iq).eq.12) then
244  ! Pas de temps adaptatif
245  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
246  if (n.GT.1) then
247  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
248  dtvr,'n=',n
249  endif
250  do indice=1,n
251  call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
252  end do
253  else if(iadv(iq).eq.13) then
254  ! Pas de temps adaptatif
255  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
256  if (n.GT.1) then
257  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
258  dtvr,'n=',n
259  endif
260  do indice=1,n
261  call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
262  end do
263  ! ----------------------------------------------------------------
264  ! Schema de pente SLOPES
265  ! ----------------------------------------------------------------
266  else if (iadv(iq).eq.20) then
267  call pentes_ini(q(1,1,iq),wg,massem,pbarug,pbarvg,0)
268 
269  ! ----------------------------------------------------------------
270  ! Schema de Prather
271  ! ----------------------------------------------------------------
272  else if (iadv(iq).eq.30) then
273  ! Pas de temps adaptatif
274  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
275  if (n.GT.1) then
276  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
277  dtvr,'n=',n
278  endif
279  call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
280  n,dtbon)
281 
282  ! ----------------------------------------------------------------
283  ! Schemas PPM Lin et Rood
284  ! ----------------------------------------------------------------
285  else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
286  iadv(iq).LE.18)) then
287 
288  ! Test sur le flux horizontal
289  ! Pas de temps adaptatif
290  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
291  if (n.GT.1) then
292  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
293  dtvr,'n=',n
294  endif
295  ! Test sur le flux vertical
296  cflmaxz=0.
297  do l=2,llm
298  do ij=iip2,ip1jm
299  aaa=wg(ij,l)*dtvr/massem(ij,l)
300  cflmaxz=max(cflmaxz,aaa)
301  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
302  cflmaxz=max(cflmaxz,bbb)
303  enddo
304  enddo
305  if (cflmaxz.GE.1) then
306  write(*,*) 'WARNING vertical','CFLmaxz=', cflmaxz
307  endif
308 
309  !-----------------------------------------------------------
310  ! Ss-prg interface LMDZ.4->PPM3d
311  !-----------------------------------------------------------
312 
313  call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
314  apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
315  unatppm,vnatppm,psppm)
316 
317  do indice=1,n
318  !----------------------------------------------------------------
319  ! VL (version PPM) horiz. et PPM vert.
320  !----------------------------------------------------------------
321  if (iadv(iq).eq.11) then
322  ! Ss-prg PPM3d de Lin
323  call ppm3d(1,qppm(1,1,iq), &
324  psppm,psppm, &
325  unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
326  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
327  fill,dum,220.)
328 
329  !-------------------------------------------------------------
330  ! Monotonic PPM
331  !-------------------------------------------------------------
332  else if (iadv(iq).eq.16) then
333  ! Ss-prg PPM3d de Lin
334  call ppm3d(1,qppm(1,1,iq), &
335  psppm,psppm, &
336  unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
337  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
338  fill,dum,220.)
339  !-------------------------------------------------------------
340 
341  !-------------------------------------------------------------
342  ! Semi Monotonic PPM
343  !-------------------------------------------------------------
344  else if (iadv(iq).eq.17) then
345  ! Ss-prg PPM3d de Lin
346  call ppm3d(1,qppm(1,1,iq), &
347  psppm,psppm, &
348  unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
349  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
350  fill,dum,220.)
351  !-------------------------------------------------------------
352 
353  !-------------------------------------------------------------
354  ! Positive Definite PPM
355  !-------------------------------------------------------------
356  else if (iadv(iq).eq.18) then
357  ! Ss-prg PPM3d de Lin
358  call ppm3d(1,qppm(1,1,iq), &
359  psppm,psppm, &
360  unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
361  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
362  fill,dum,220.)
363  !-------------------------------------------------------------
364  endif
365  enddo
366  !-----------------------------------------------------------------
367  ! Ss-prg interface PPM3d-LMDZ.4
368  !-----------------------------------------------------------------
369  call interpost(q(1,1,iq),qppm(1,1,iq))
370  endif
371  !----------------------------------------------------------------------
372 
373  !-----------------------------------------------------------------
374  ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
375  ! et Nord j=1
376  !-----------------------------------------------------------------
377 
378  ! call traceurpole(q(1,1,iq),massem)
379 
380  ! calcul du temps cpu pour un schema donne
381 
382  ! call clock(t_final)
383  !ym tps_cpu=t_final-t_initial
384  !ym cpuadv(iq)=cpuadv(iq)+tps_cpu
385 
386  end DO
387 
388 
389  !------------------------------------------------------------------
390  ! on reinitialise a zero les flux de masse cumules
391  !---------------------------------------------------
392  iadvtr=0
393 
394  ENDIF ! if iadvtr.EQ.iapp_tracvl
395 
396 END SUBROUTINE advtrac