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

Last change on this file since 965 was 960, checked in by lsce, 16 years ago
  • Ajoute du parametre config_inca dans conf_gcm.F config_inca='none'(sans INCA, par defaut) config_inca='chem'(avec INCA config chemie) config_inca='aero'(avec INCA config aerosol)
  • Menage parmis les cles CPP INCA
  • Enleve le calcul d'omega dans calfis.F et active le calcul correspondant dans advtrac.F(avant uniquement pour INCA).

JG

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 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 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
163c$OMP BARRIER
164        CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
165c$OMP BARRIER
166
167c$OMP BARRIER
168c$OMP MASTER     
169      p_tmp(ijb:ije,1:llmp1)=p(ijb:ije,1:llmp1)
170      pk_tmp(ijb:ije,1:llm)=pk(ijb:ije,1:llm)
171      teta_tmp(ijb:ije,1:llm)=teta(ijb:ije,1:llm)
172      call VTb(VTHallo)
173      call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
174     *                             jj_Nb_vanleer,0,0,Request_vanleer)
175      call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
176     *                             jj_Nb_vanleer,1,0,Request_vanleer)
177      call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
178     *                             jj_Nb_vanleer,0,0,Request_vanleer)
179      call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
180     *                             jj_Nb_vanleer,0,0,Request_vanleer)
181      call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
182     *                             jj_Nb_vanleer,1,1,Request_vanleer)
183      call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
184     *                             jj_Nb_vanleer,1,1,Request_vanleer)
185      call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
186     *                             jj_Nb_vanleer,1,1,Request_vanleer)
187      do j=1,nqmx
188        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
189     *                             jj_nb_vanleer,0,0,Request_vanleer)
190      enddo
191     
192      call SetDistrib(jj_nb_vanleer)
193      call SendRequest(Request_vanleer)
194      call WaitRequest(Request_vanleer)
195
196      call VTe(VTHallo)
197     
198      call VTb(VTadvection)
199      call start_timer(timer_vanleer)
200c$OMP END MASTER
201c$OMP BARRIER
202     
203      ! ... Flux de masse diaganostiques traceurs
204         ijb=ij_begin
205         ije=ij_end
206         flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/FLOAT(iapp_tracvl)
207
208c  test sur l'eventuelle creation de valeurs negatives de la masse
209         ijb=ij_begin
210         ije=ij_end
211         if (pole_nord) ijb=ij_begin+iip1
212         if (pole_sud) ije=ij_end-iip1
213         
214c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
215         DO l=1,llm-1
216            DO ij = ijb+1,ije
217              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
218     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
219     s                  +       wg(ij,l+1)  - wg(ij,l)
220            ENDDO
221           
222c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
223c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
224           
225            do ij=ijb,ije-iip1+1,iip1
226              zdp(ij)=zdp(ij+iip1-1)
227            enddo
228           
229            DO ij = ijb,ije
230               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
231            ENDDO
232
233
234c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
235c  ym ---> eventuellement a revoir
236            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
237           
238            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
239            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
240     s        '   MAX:', zdpmax
241            ENDIF
242
243         ENDDO
244c$OMP END DO NOWAIT
245
246c-------------------------------------------------------------------
247c   Advection proprement dite (Modification Le Croller (07/2001)
248c-------------------------------------------------------------------
249
250c----------------------------------------------------
251c        Calcul des moyennes basées sur la masse
252c----------------------------------------------------
253
254cym      ----> Normalement, inutile pour les schémas classiques
255cym      ----> Revérifier lors de la parallélisation des autres schemas
256   
257cym          call massbar_p(massem,massebx,masseby)         
258
259          call vlspltgen_p( q,iadv, 2., massem, wg ,
260     *                    pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
261
262         
263          GOTO 1234     
264c-----------------------------------------------------------
265c     Appel des sous programmes d'advection
266c-----------------------------------------------------------
267      do iq=1,nqmx
268c        call clock(t_initial)
269        if(iadv(iq) == 0) cycle
270c   ----------------------------------------------------------------
271c   Schema de Van Leer I MUSCL
272c   ----------------------------------------------------------------
273        if(iadv(iq).eq.10) THEN
274     
275            call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
276
277c   ----------------------------------------------------------------
278c   Schema "pseudo amont" + test sur humidite specifique
279C    pour la vapeur d'eau. F. Codron
280c   ----------------------------------------------------------------
281        else if(iadv(iq).eq.14) then
282c
283cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
284           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
285     *                 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
286c   ----------------------------------------------------------------
287c   Schema de Frederic Hourdin
288c   ----------------------------------------------------------------
289        else if(iadv(iq).eq.12) then
290          stop 'advtrac : schema non parallelise'
291c            Pas de temps adaptatif
292           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
293           if (n.GT.1) then
294           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
295     s             dtvr,'n=',n
296           endif
297           do indice=1,n
298            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
299           end do
300        else if(iadv(iq).eq.13) 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,2)
310          end do
311c   ----------------------------------------------------------------
312c   Schema de pente SLOPES
313c   ----------------------------------------------------------------
314        else if (iadv(iq).eq.20) then
315          stop 'advtrac : schema non parallelise'
316
317            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
318
319c   ----------------------------------------------------------------
320c   Schema de Prather
321c   ----------------------------------------------------------------
322        else if (iadv(iq).eq.30) then
323          stop 'advtrac : schema non parallelise'
324c            Pas de temps adaptatif
325           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
326           if (n.GT.1) then
327           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
328     s             dtvr,'n=',n
329           endif
330           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
331     s                     n,dtbon)
332c   ----------------------------------------------------------------
333c   Schemas PPM Lin et Rood
334c   ----------------------------------------------------------------
335       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
336     s                     iadv(iq).LE.18)) then
337
338           stop 'advtrac : schema non parallelise'
339
340c        Test sur le flux horizontal
341c        Pas de temps adaptatif
342         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
343         if (n.GT.1) then
344           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
345     s             dtvr,'n=',n
346         endif
347c        Test sur le flux vertical
348         CFLmaxz=0.
349         do l=2,llm
350           do ij=iip2,ip1jm
351            aaa=wg(ij,l)*dtvr/massem(ij,l)
352            CFLmaxz=max(CFLmaxz,aaa)
353            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
354            CFLmaxz=max(CFLmaxz,bbb)
355           enddo
356         enddo
357         if (CFLmaxz.GE.1) then
358            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
359         endif
360
361c-----------------------------------------------------------
362c        Ss-prg interface LMDZ.4->PPM3d
363c-----------------------------------------------------------
364
365          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
366     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
367     s                 unatppm,vnatppm,psppm)
368
369          do indice=1,n
370c---------------------------------------------------------------------
371c                         VL (version PPM) horiz. et PPM vert.
372c---------------------------------------------------------------------
373                if (iadv(iq).eq.11) then
374c                  Ss-prg PPM3d de Lin
375                  call ppm3d(1,qppm(1,1,iq),
376     s                       psppm,psppm,
377     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
378     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
379     s                       fill,dum,220.)
380
381c----------------------------------------------------------------------
382c                           Monotonic PPM
383c----------------------------------------------------------------------
384               else if (iadv(iq).eq.16) then
385c                  Ss-prg PPM3d de Lin
386                  call ppm3d(1,qppm(1,1,iq),
387     s                       psppm,psppm,
388     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
389     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
390     s                       fill,dum,220.)
391c---------------------------------------------------------------------
392
393c---------------------------------------------------------------------
394c                           Semi Monotonic PPM
395c---------------------------------------------------------------------
396               else if (iadv(iq).eq.17) then
397c                  Ss-prg PPM3d de Lin
398                  call ppm3d(1,qppm(1,1,iq),
399     s                       psppm,psppm,
400     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
401     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
402     s                       fill,dum,220.)
403c---------------------------------------------------------------------
404
405c---------------------------------------------------------------------
406c                         Positive Definite PPM
407c---------------------------------------------------------------------
408                else if (iadv(iq).eq.18) then
409c                  Ss-prg PPM3d de Lin
410                  call ppm3d(1,qppm(1,1,iq),
411     s                       psppm,psppm,
412     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
413     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
414     s                       fill,dum,220.)
415c---------------------------------------------------------------------
416                endif
417            enddo
418c-----------------------------------------------------------------
419c               Ss-prg interface PPM3d-LMDZ.4
420c-----------------------------------------------------------------
421                  call interpost(q(1,1,iq),qppm(1,1,iq))
422            endif
423c----------------------------------------------------------------------
424
425c-----------------------------------------------------------------
426c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
427c et Nord j=1
428c-----------------------------------------------------------------
429
430c                  call traceurpole(q(1,1,iq),massem)
431
432c calcul du temps cpu pour un schema donne
433
434c                  call clock(t_final)
435cym                  tps_cpu=t_final-t_initial
436cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
437
438       end DO
439
4401234  CONTINUE
441c$OMP BARRIER
442c$OMP MASTER
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
453       CALL qminimum_p( q, 2, finmasse )
454
455c------------------------------------------------------------------
456c   on reinitialise a zero les flux de masse cumules
457c---------------------------------------------------
458c          iadvtr=0
459        call VTe(VTadvection)
460        call stop_timer(timer_vanleer)
461
462        call VTb(VThallo)
463        do j=1,nqmx
464          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
465     *                             jj_nb_caldyn,0,0,Request_vanleer)
466        enddo
467
468        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
469     *       jj_nb_caldyn,0,0,Request_vanleer)
470
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)
477c$OMP END MASTER
478c$OMP BARRIER   
479          iadvtr=0
480       ENDIF ! if iadvtr.EQ.iapp_tracvl
481
482       RETURN
483       END
484
Note: See TracBrowser for help on using the repository browser.