My Project
 All Classes Files Functions Variables Macros
advtrac_p.F90
Go to the documentation of this file.
1 ! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $
2 
3 SUBROUTINE advtrac_p(pbaru,pbarv , p, masse,q,iapptrac,teta, flxw, pk)
4 
5  ! Auteur : F. Hourdin
6  !
7  ! Modif. P. Le Van (20/12/97)
8  ! F. Codron (10/99)
9  ! D. Le Croller (07/2001)
10  ! M.A Filiberti (04/2002)
11  !
12  USE parallel
13  USE write_field_p
14  USE bands
15  USE mod_hallo
16  USE vampir
17  USE times
18  USE infotrac
19  USE control_mod
20  IMPLICIT NONE
21  !
22  include "dimensions.h"
23  include "paramet.h"
24  include "comconst.h"
25  include "comvert.h"
26  include "comdissip.h"
27  include "comgeom2.h"
28  include "logic.h"
29  include "temps.h"
30  include "ener.h"
31  include "description.h"
32 
33  !-------------------------------------------------------------------
34  ! Arguments
35  !-------------------------------------------------------------------
36  ! Ajout PPM
37  !--------------------------------------------------------
38  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
39  !--------------------------------------------------------
40  INTEGER iapptrac
41  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
42  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
43  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
44  REAL pk(ip1jmp1,llm)
45  REAL :: flxw(ip1jmp1,llm)
46 
47  !-------------------------------------------------------------
48  ! Variables locales
49  !-------------------------------------------------------------
50 
51  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
52  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
53  REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
54  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
55  INTEGER iadvtr
56  INTEGER ij,l,iq,iiq
57  REAL zdpmin, zdpmax
58  SAVE iadvtr, massem, pbaruc, pbarvc
59  DATA iadvtr/0/
60  !$OMP THREADPRIVATE(iadvtr)
61  !----------------------------------------------------------
62  ! Rajouts pour PPM
63  !----------------------------------------------------------
64  INTEGER indice,n
65  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
66  REAL cflmaxz,aaa,bbb ! CFL maximum
67  REAL psppm(iim,jjp1) ! pression au sol
68  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
69  REAL qppm(iim*jjp1,llm,nqtot)
70  REAL fluxwppm(iim,jjp1,llm)
71  REAL apppm(llmp1), bpppm(llmp1)
72  LOGICAL dum,fill
73  DATA fill/.true./
74  DATA dum/.true./
75  REAL,SAVE :: finmasse(ip1jmp1,llm)
76  integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
77  type(request) :: request_vanleer
78  REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
79  REAL,SAVE :: teta_tmp(ip1jmp1,llm)
80  REAL,SAVE :: pk_tmp(ip1jmp1,llm)
81 
82  ijb_u=ij_begin
83  ije_u=ij_end
84 
85  ijb_v=ij_begin-iip1
86  ije_v=ij_end
87  if (pole_nord) ijb_v=ij_begin
88  if (pole_sud) ije_v=ij_end-iip1
89 
90  IF(iadvtr.EQ.0) THEN
91  ! CALL initial0(ijp1llm,pbaruc)
92  ! CALL initial0(ijmllm,pbarvc)
93  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
94  DO l=1,llm
95  pbaruc(ijb_u:ije_u,l)=0.
96  pbarvc(ijb_v:ije_v,l)=0.
97  ENDDO
98  !$OMP END DO NOWAIT
99  ENDIF
100 
101  ! accumulation des flux de masse horizontaux
102  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
103  DO l=1,llm
104  DO ij = ijb_u,ije_u
105  pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
106  ENDDO
107  DO ij = ijb_v,ije_v
108  pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
109  ENDDO
110  ENDDO
111  !$OMP END DO NOWAIT
112 
113  ! selection de la masse instantannee des mailles avant le transport.
114  IF(iadvtr.EQ.0) THEN
115 
116  ! CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
117  ijb=ij_begin
118  ije=ij_end
119 
120  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
121  DO l=1,llm
122  massem(ijb:ije,l)=masse(ijb:ije,l)
123  ENDDO
124  !$OMP END DO NOWAIT
125 
126  !cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
127  !
128  ENDIF ! of IF(iadvtr.EQ.0)
129 
130  iadvtr = iadvtr+1
131 
132  !$OMP MASTER
133  iapptrac = iadvtr
134  !$OMP END MASTER
135 
136  ! Test pour savoir si on advecte a ce pas de temps
137 
138  IF ( iadvtr.EQ.iapp_tracvl ) THEN
139  !$OMP MASTER
140  call suspend_timer(timer_caldyn)
141  !$OMP END MASTER
142 
143  ijb=ij_begin
144  ije=ij_end
145 
146 
147  !c .. Modif P.Le Van ( 20/12/97 ) ....
148  !c
149 
150  ! traitement des flux de masse avant advection.
151  ! 1. calcul de w
152  ! 2. groupement des mailles pres du pole.
153 
154  CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
155 
156  !$OMP BARRIER
157 
158  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
159  DO l=1,llmp1
160  p_tmp(ijb:ije,l)=p(ijb:ije,l)
161  ENDDO
162  !$OMP END DO NOWAIT
163 
164  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
165  DO l=1,llm
166  pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
167  teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
168  ENDDO
169  !$OMP END DO NOWAIT
170 
171  !$OMP MASTER
172  call vtb(vthallo)
173  !$OMP END MASTER
174 
175  call register_swapfieldhallo(pbarug,pbarug,ip1jmp1,llm, &
176  jj_nb_vanleer,0,0,request_vanleer)
177  call register_swapfieldhallo(pbarvg,pbarvg,ip1jm,llm, &
178  jj_nb_vanleer,1,0,request_vanleer)
179  call register_swapfieldhallo(massem,massem,ip1jmp1,llm, &
180  jj_nb_vanleer,0,0,request_vanleer)
181  call register_swapfieldhallo(wg,wg,ip1jmp1,llm, &
182  jj_nb_vanleer,0,0,request_vanleer)
183  call register_swapfieldhallo(teta_tmp,teta_tmp,ip1jmp1,llm, &
184  jj_nb_vanleer,1,1,request_vanleer)
185  call register_swapfieldhallo(p_tmp,p_tmp,ip1jmp1,llmp1, &
186  jj_nb_vanleer,1,1,request_vanleer)
187  call register_swapfieldhallo(pk_tmp,pk_tmp,ip1jmp1,llm, &
188  jj_nb_vanleer,1,1,request_vanleer)
189  do j=1,nqtot
190  call register_swapfieldhallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
191  jj_nb_vanleer,0,0,request_vanleer)
192  enddo
193 
194  call sendrequest(request_vanleer)
195  !$OMP BARRIER
196  call waitrequest(request_vanleer)
197 
198 
199  !$OMP BARRIER
200  !$OMP MASTER
201  call setdistrib(jj_nb_vanleer)
202  call vte(vthallo)
203  call vtb(vtadvection)
204  call start_timer(timer_vanleer)
205  !$OMP END MASTER
206  !$OMP BARRIER
207 
208  ! ... Flux de masse diaganostiques traceurs
209  ijb=ij_begin
210  ije=ij_end
211  flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
212 
213  ! test sur l'eventuelle creation de valeurs negatives de la masse
214  ijb=ij_begin
215  ije=ij_end
216  if (pole_nord) ijb=ij_begin+iip1
217  if (pole_sud) ije=ij_end-iip1
218 
219  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
220  DO l=1,llm-1
221  DO ij = ijb+1,ije
222  zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
223  - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
224  + wg(ij,l+1) - wg(ij,l)
225  ENDDO
226 
227  ! CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
228  ! ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
229 
230  do ij=ijb,ije-iip1+1,iip1
231  zdp(ij)=zdp(ij+iip1-1)
232  enddo
233 
234  DO ij = ijb,ije
235  zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
236  ENDDO
237 
238 
239  ! CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
240  ! ym ---> eventuellement a revoir
241  CALL minmax( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
242 
243  IF(max(abs(zdpmin),abs(zdpmax)).GT.0.5) THEN
244  print*,'WARNING DP/P l=',l,' MIN:',zdpmin, &
245  ' MAX:', zdpmax
246  ENDIF
247 
248  ENDDO
249  !$OMP END DO NOWAIT
250 
251  !-------------------------------------------------------------------
252  ! Advection proprement dite (Modification Le Croller (07/2001)
253  !-------------------------------------------------------------------
254 
255  !----------------------------------------------------
256  ! Calcul des moyennes basées sur la masse
257  !----------------------------------------------------
258 
259  !ym ----> Normalement, inutile pour les schémas classiques
260  !ym ----> Revérifier lors de la parallélisation des autres schemas
261 
262  !ym call massbar_p(massem,massebx,masseby)
263 
264  call vlspltgen_p( q,iadv, 2., massem, wg , &
265  pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
266 
267 
268  goto 1234
269  !-----------------------------------------------------------
270  ! Appel des sous programmes d'advection
271  !-----------------------------------------------------------
272  do iq=1,nqtot
273  ! call clock(t_initial)
274  if(iadv(iq) == 0) cycle
275  ! ----------------------------------------------------------------
276  ! Schema de Van Leer I MUSCL
277  ! ----------------------------------------------------------------
278  if(iadv(iq).eq.10) THEN
279 
280  call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
281 
282  ! ----------------------------------------------------------------
283  ! Schema "pseudo amont" + test sur humidite specifique
284  ! pour la vapeur d'eau. F. Codron
285  ! ----------------------------------------------------------------
286  else if(iadv(iq).eq.14) then
287  !
288  !ym stop 'advtrac : appel à vlspltqs :schema non parallelise'
289  CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
290  pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
291  ! ----------------------------------------------------------------
292  ! Schema de Frederic Hourdin
293  ! ----------------------------------------------------------------
294  else if(iadv(iq).eq.12) then
295  stop 'advtrac : schema non parallelise'
296  ! Pas de temps adaptatif
297  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
298  if (n.GT.1) then
299  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
300  dtvr,'n=',n
301  endif
302  do indice=1,n
303  call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
304  end do
305  else if(iadv(iq).eq.13) then
306  stop 'advtrac : schema non parallelise'
307  ! Pas de temps adaptatif
308  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
309  if (n.GT.1) then
310  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
311  dtvr,'n=',n
312  endif
313  do indice=1,n
314  call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
315  end do
316  ! ----------------------------------------------------------------
317  ! Schema de pente SLOPES
318  ! ----------------------------------------------------------------
319  else if (iadv(iq).eq.20) then
320  stop 'advtrac : schema non parallelise'
321 
322  call pentes_ini(q(1,1,iq),wg,massem,pbarug,pbarvg,0)
323 
324  ! ----------------------------------------------------------------
325  ! Schema de Prather
326  ! ----------------------------------------------------------------
327  else if (iadv(iq).eq.30) then
328  stop 'advtrac : schema non parallelise'
329  ! Pas de temps adaptatif
330  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
331  if (n.GT.1) then
332  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
333  dtvr,'n=',n
334  endif
335  call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
336  n,dtbon)
337  ! ----------------------------------------------------------------
338  ! Schemas PPM Lin et Rood
339  ! ----------------------------------------------------------------
340  else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
341  iadv(iq).LE.18)) then
342 
343  stop 'advtrac : schema non parallelise'
344 
345  ! Test sur le flux horizontal
346  ! Pas de temps adaptatif
347  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
348  if (n.GT.1) then
349  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
350  dtvr,'n=',n
351  endif
352  ! Test sur le flux vertical
353  cflmaxz=0.
354  do l=2,llm
355  do ij=iip2,ip1jm
356  aaa=wg(ij,l)*dtvr/massem(ij,l)
357  cflmaxz=max(cflmaxz,aaa)
358  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
359  cflmaxz=max(cflmaxz,bbb)
360  enddo
361  enddo
362  if (cflmaxz.GE.1) then
363  write(*,*) 'WARNING vertical','CFLmaxz=', cflmaxz
364  endif
365 
366  !-----------------------------------------------------------
367  ! Ss-prg interface LMDZ.4->PPM3d
368  !-----------------------------------------------------------
369 
370  call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
371  apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
372  unatppm,vnatppm,psppm)
373 
374  do indice=1,n
375  !----------------------------------------------------------------
376  ! VL (version PPM) horiz. et PPM vert.
377  !----------------------------------------------------------------
378  if (iadv(iq).eq.11) then
379  ! Ss-prg PPM3d de Lin
380  call ppm3d(1,qppm(1,1,iq), &
381  psppm,psppm, &
382  unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
383  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
384  fill,dum,220.)
385 
386  !-------------------------------------------------------------
387  ! Monotonic PPM
388  !-------------------------------------------------------------
389  else if (iadv(iq).eq.16) then
390  ! Ss-prg PPM3d de Lin
391  call ppm3d(1,qppm(1,1,iq), &
392  psppm,psppm, &
393  unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
394  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
395  fill,dum,220.)
396  !-------------------------------------------------------------
397 
398  !-------------------------------------------------------------
399  ! Semi Monotonic PPM
400  !-------------------------------------------------------------
401  else if (iadv(iq).eq.17) then
402  ! Ss-prg PPM3d de Lin
403  call ppm3d(1,qppm(1,1,iq), &
404  psppm,psppm, &
405  unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
406  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
407  fill,dum,220.)
408  !-------------------------------------------------------------
409 
410  !-------------------------------------------------------------
411  ! Positive Definite PPM
412  !-------------------------------------------------------------
413  else if (iadv(iq).eq.18) then
414  ! Ss-prg PPM3d de Lin
415  call ppm3d(1,qppm(1,1,iq), &
416  psppm,psppm, &
417  unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
418  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
419  fill,dum,220.)
420  !-------------------------------------------------------------
421  endif
422  enddo
423  !-----------------------------------------------------------------
424  ! Ss-prg interface PPM3d-LMDZ.4
425  !-----------------------------------------------------------------
426  call interpost(q(1,1,iq),qppm(1,1,iq))
427  endif
428  !----------------------------------------------------------------------
429 
430  !-----------------------------------------------------------------
431  ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
432  ! et Nord j=1
433  !-----------------------------------------------------------------
434 
435  ! call traceurpole(q(1,1,iq),massem)
436 
437  ! calcul du temps cpu pour un schema donne
438 
439  ! call clock(t_final)
440  !ym tps_cpu=t_final-t_initial
441  !ym cpuadv(iq)=cpuadv(iq)+tps_cpu
442 
443  end DO
444 
445 1234 CONTINUE
446  !$OMP BARRIER
447 
448  if (planet_type=="earth") then
449 
450  ijb=ij_begin
451  ije=ij_end
452 
453  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
454  DO l = 1, llm
455  DO ij = ijb, ije
456  finmasse(ij,l) = p(ij,l) - p(ij,l+1)
457  ENDDO
458  ENDDO
459  !$OMP END DO
460 
461  CALL qminimum_p( q, 2, finmasse )
462 
463  !------------------------------------------------------------------
464  ! on reinitialise a zero les flux de masse cumules
465  !---------------------------------------------------
466  ! iadvtr=0
467 
468  !$OMP MASTER
469  call vte(vtadvection)
470  call stop_timer(timer_vanleer)
471  call vtb(vthallo)
472  !$OMP END MASTER
473 
474  do j=1,nqtot
475  call register_swapfieldhallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
476  jj_nb_caldyn,0,0,request_vanleer)
477  enddo
478 
479  call register_swapfieldhallo(flxw,flxw,ip1jmp1,llm, &
480  jj_nb_caldyn,0,0,request_vanleer)
481 
482  call sendrequest(request_vanleer)
483  !$OMP BARRIER
484  call waitrequest(request_vanleer)
485 
486  !$OMP BARRIER
487  !$OMP MASTER
488  call setdistrib(jj_nb_caldyn)
489  call vte(vthallo)
490  call resume_timer(timer_caldyn)
491  !$OMP END MASTER
492  !$OMP BARRIER
493  iadvtr=0
494  endif ! of if (planet_type=="earth")
495  ENDIF ! if iadvtr.EQ.iapp_tracvl
496 
497 END SUBROUTINE advtrac_p