source: LMDZ4/branches/V3_test/libf/dyn3dpar/advtrac_p.F @ 2286

Last change on this file since 2286 was 727, checked in by Laurent Fairhead, 18 years ago

Modif pour rendre KLDZ et INCA plus indépendants ACo
LF

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