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

Last change on this file since 654 was 630, checked in by Laurent Fairhead, 20 years ago

Import d'une version parallele de la dynamique YM
LF

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