source: trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90 @ 1000

Last change on this file since 1000 was 953, checked in by slebonnois, 12 years ago

SL: optimisation pour le parallèle suite à tests Venus / petite correction appels routines secondaires dans Venus et Titan

File size: 17.0 KB
Line 
1! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $
2
3SUBROUTINE 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 instantanee 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 diagnostiques 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!!! ATTENTION !!!! TOUT CE QUI EST ENTRE ICI ET 1234 EST OBSOLETE !!!!!!!
269     GOTO 1234
270!!! ATTENTION !!!!
271
272     !-----------------------------------------------------------
273     !     Appel des sous programmes d'advection
274     !-----------------------------------------------------------
275     do iq=1,nqtot
276        !        call clock(t_initial)
277        if(iadv(iq) == 0) cycle
278        !   ----------------------------------------------------------------
279        !   Schema de Van Leer I MUSCL
280        !   ----------------------------------------------------------------
281        if(iadv(iq).eq.10) THEN
282
283           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
284
285           !   ----------------------------------------------------------------
286           !   Schema "pseudo amont" + test sur humidite specifique
287           !    pour la vapeur d'eau. F. Codron
288           !   ----------------------------------------------------------------
289        else if(iadv(iq).eq.14) then
290           !
291           !ym        stop 'advtrac : appel à vlspltqs :schema non parallelise'
292           CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
293                pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
294           !   ----------------------------------------------------------------
295           !   Schema de Frederic Hourdin
296           !   ----------------------------------------------------------------
297        else if(iadv(iq).eq.12) then
298           stop 'advtrac : schema non parallelise'
299           !            Pas de temps adaptatif
300           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
301           if (n.GT.1) then
302              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
303                   dtvr,'n=',n
304           endif
305           do indice=1,n
306              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
307           end do
308        else if(iadv(iq).eq.13) then
309           stop 'advtrac : schema non parallelise'
310           !            Pas de temps adaptatif
311           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
312           if (n.GT.1) then
313              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
314                   dtvr,'n=',n
315           endif
316           do indice=1,n
317              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
318           end do
319           !   ----------------------------------------------------------------
320           !   Schema de pente SLOPES
321           !   ----------------------------------------------------------------
322        else if (iadv(iq).eq.20) then
323           stop 'advtrac : schema non parallelise'
324
325           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
326
327           !   ----------------------------------------------------------------
328           !   Schema de Prather
329           !   ----------------------------------------------------------------
330        else if (iadv(iq).eq.30) then
331           stop 'advtrac : schema non parallelise'
332           !            Pas de temps adaptatif
333           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
334           if (n.GT.1) then
335              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
336                   dtvr,'n=',n
337           endif
338           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
339                n,dtbon)
340           !   ----------------------------------------------------------------
341           !   Schemas PPM Lin et Rood
342           !   ----------------------------------------------------------------
343        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
344             iadv(iq).LE.18)) then
345
346           stop 'advtrac : schema non parallelise'
347
348           !        Test sur le flux horizontal
349           !        Pas de temps adaptatif
350           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
351           if (n.GT.1) then
352              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
353                   dtvr,'n=',n
354           endif
355           !        Test sur le flux vertical
356           CFLmaxz=0.
357           do l=2,llm
358              do ij=iip2,ip1jm
359                 aaa=wg(ij,l)*dtvr/massem(ij,l)
360                 CFLmaxz=max(CFLmaxz,aaa)
361                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
362                 CFLmaxz=max(CFLmaxz,bbb)
363              enddo
364           enddo
365           if (CFLmaxz.GE.1) then
366              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
367           endif
368
369           !-----------------------------------------------------------
370           !        Ss-prg interface LMDZ.4->PPM3d
371           !-----------------------------------------------------------
372
373           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
374                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
375                unatppm,vnatppm,psppm)
376
377           do indice=1,n
378              !----------------------------------------------------------------
379              !                         VL (version PPM) horiz. et PPM vert.
380              !----------------------------------------------------------------
381              if (iadv(iq).eq.11) then
382                 !                  Ss-prg PPM3d de Lin
383                 call ppm3d(1,qppm(1,1,iq), &
384                      psppm,psppm, &
385                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
386                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
387                      fill,dum,220.)
388
389                 !-------------------------------------------------------------
390                 !                           Monotonic PPM
391                 !-------------------------------------------------------------
392              else if (iadv(iq).eq.16) then
393                 !                  Ss-prg PPM3d de Lin
394                 call ppm3d(1,qppm(1,1,iq), &
395                      psppm,psppm, &
396                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
397                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
398                      fill,dum,220.)
399                 !-------------------------------------------------------------
400
401                 !-------------------------------------------------------------
402                 !                           Semi Monotonic PPM
403                 !-------------------------------------------------------------
404              else if (iadv(iq).eq.17) then
405                 !                  Ss-prg PPM3d de Lin
406                 call ppm3d(1,qppm(1,1,iq), &
407                      psppm,psppm, &
408                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
409                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
410                      fill,dum,220.)
411                 !-------------------------------------------------------------
412
413                 !-------------------------------------------------------------
414                 !                         Positive Definite PPM
415                 !-------------------------------------------------------------
416              else if (iadv(iq).eq.18) then
417                 !                  Ss-prg PPM3d de Lin
418                 call ppm3d(1,qppm(1,1,iq), &
419                      psppm,psppm, &
420                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
421                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
422                      fill,dum,220.)
423                 !-------------------------------------------------------------
424              endif
425           enddo
426           !-----------------------------------------------------------------
427           !               Ss-prg interface PPM3d-LMDZ.4
428           !-----------------------------------------------------------------
429           call interpost(q(1,1,iq),qppm(1,1,iq))
430        endif
431        !----------------------------------------------------------------------
432
433        !-----------------------------------------------------------------
434        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
435        ! et Nord j=1
436        !-----------------------------------------------------------------
437
438        !                  call traceurpole(q(1,1,iq),massem)
439
440        ! calcul du temps cpu pour un schema donne
441
442        !                  call clock(t_final)
443        !ym                  tps_cpu=t_final-t_initial
444        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
445
446     end DO
447
448!!! ATTENTION !!!!
4491234 CONTINUE
450!!! ATTENTION !!!! LE CODE REPREND ICI !!!!!!!!
451
452     !$OMP BARRIER
453
454     if (planet_type=="earth") then
455
456        ijb=ij_begin
457        ije=ij_end
458
459        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
460        DO l = 1, llm
461           DO ij = ijb, ije
462              finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
463           ENDDO
464        ENDDO
465        !$OMP END DO
466
467        CALL qminimum_p( q, 2, finmasse )
468
469     endif ! of if (planet_type=="earth")
470
471        !------------------------------------------------------------------
472        !   on reinitialise a zero les flux de masse cumules
473        !---------------------------------------------------
474        !          iadvtr=0
475
476        !$OMP MASTER
477        call VTe(VTadvection)
478        call stop_timer(timer_vanleer)
479        call VTb(VThallo)
480        !$OMP END MASTER
481
482
483        do j=1,nqtot
484           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
485                jj_nb_caldyn,0,0,Request_vanleer)
486        enddo
487
488        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
489             jj_nb_caldyn,0,0,Request_vanleer)
490
491        call SendRequest(Request_vanleer)
492        !$OMP BARRIER
493        call WaitRequest(Request_vanleer)     
494
495        !$OMP BARRIER
496        !$OMP MASTER
497        call SetDistrib(jj_nb_caldyn)
498        call VTe(VThallo)
499        call resume_timer(timer_caldyn)
500 !$OMP END MASTER
501 !$OMP BARRIER 
502        iadvtr=0
503  ENDIF ! if iadvtr.EQ.iapp_tracvl
504
505END SUBROUTINE advtrac_p
Note: See TracBrowser for help on using the repository browser.