source: LMDZ4/trunk/libf/dyn3dpar/advtrac_p.F @ 1098

Last change on this file since 1098 was 985, checked in by Laurent Fairhead, 16 years ago

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF

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