source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/advtrac.F @ 1151

Last change on this file since 1151 was 1114, checked in by jghattas, 16 years ago

Creation du module infotrac:

  • contient les variables de advtrac.h
  • contient la subroutine iniadvtrac renommer en infotrac_init
  • le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
  • ce module est utilise dans la dynamique et la physique
  • contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.6 KB
RevLine 
[524]1!
2! $Header$
3!
4c
5c
6      SUBROUTINE advtrac(pbaru,pbarv ,
7     *                   p,  masse,q,iapptrac,teta,
8     *                  flxw,
[762]9     *                  pk)
[524]10c     Auteur :  F. Hourdin
11c
12c     Modif. P. Le Van     (20/12/97)
13c            F. Codron     (10/99)
14c            D. Le Croller (07/2001)
15c            M.A Filiberti (04/2002)
16c
[1114]17      USE infotrac
18
[524]19      IMPLICIT NONE
20c
21#include "dimensions.h"
22#include "paramet.h"
23#include "comconst.h"
24#include "comvert.h"
25#include "comdissip.h"
26#include "comgeom2.h"
27#include "logic.h"
28#include "temps.h"
29#include "control.h"
30#include "ener.h"
31#include "description.h"
32
33c-------------------------------------------------------------------
34c     Arguments
35c-------------------------------------------------------------------
36c     Ajout PPM
37c--------------------------------------------------------
38      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
39c--------------------------------------------------------
40      INTEGER iapptrac
41      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
[1114]42      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
[524]43      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
44      REAL pk(ip1jmp1,llm)
[960]45      REAL flxw(ip1jmp1,llm)
[524]46
47c-------------------------------------------------------------
48c     Variables locales
49c-------------------------------------------------------------
50
51      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
52      REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
53      REAL 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      EXTERNAL  minmax
59      SAVE iadvtr, massem, pbaruc, pbarvc
60      DATA iadvtr/0/
61c----------------------------------------------------------
62c     Rajouts pour PPM
63c----------------------------------------------------------
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)
[1114]69      REAL qppm(iim*jjp1,llm,nqtot)
[524]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
76
77      IF(iadvtr.EQ.0) THEN
78         CALL initial0(ijp1llm,pbaruc)
79         CALL initial0(ijmllm,pbarvc)
80      ENDIF
81
82c   accumulation des flux de masse horizontaux
83      DO l=1,llm
84         DO ij = 1,ip1jmp1
85            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
86         ENDDO
87         DO ij = 1,ip1jm
88            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
89         ENDDO
90      ENDDO
91
92c   selection de la masse instantannee des mailles avant le transport.
93      IF(iadvtr.EQ.0) THEN
94
95         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
96ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
97c
98      ENDIF
99
100      iadvtr   = iadvtr+1
101      iapptrac = iadvtr
102
103
104c   Test pour savoir si on advecte a ce pas de temps
105      IF ( iadvtr.EQ.iapp_tracvl ) THEN
106
107cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
108cc
109
110c   traitement des flux de masse avant advection.
111c     1. calcul de w
112c     2. groupement des mailles pres du pole.
113
114        CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
115
116      ! ... Flux de masse diaganostiques traceurs
117      flxw = wg / FLOAT(iapp_tracvl)
118
119c  test sur l'eventuelle creation de valeurs negatives de la masse
120         DO l=1,llm-1
121            DO ij = iip2+1,ip1jm
122              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
123     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
124     s                  +       wg(ij,l+1)  - wg(ij,l)
125            ENDDO
126            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
127            DO ij = iip2,ip1jm
128               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
129            ENDDO
130
131
132            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
133
134            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
135            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
136     s        '   MAX:', zdpmax
137            ENDIF
138
139         ENDDO
140
141c-------------------------------------------------------------------
142c   Advection proprement dite (Modification Le Croller (07/2001)
143c-------------------------------------------------------------------
144
145c----------------------------------------------------
146c        Calcul des moyennes basées sur la masse
147c----------------------------------------------------
148          call massbar(massem,massebx,masseby)         
149
150c-----------------------------------------------------------
151c     Appel des sous programmes d'advection
152c-----------------------------------------------------------
[1114]153      do iq=1,nqtot
[524]154c        call clock(t_initial)
155        if(iadv(iq) == 0) cycle
156c   ----------------------------------------------------------------
157c   Schema de Van Leer I MUSCL
158c   ----------------------------------------------------------------
159        if(iadv(iq).eq.10) THEN
160            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
161
162
163c   ----------------------------------------------------------------
164c   Schema "pseudo amont" + test sur humidite specifique
165C    pour la vapeur d'eau. F. Codron
166c   ----------------------------------------------------------------
167        else if(iadv(iq).eq.14) then
168c
169           CALL vlspltqs( q(1,1,1), 2., massem, wg ,
170     *                 pbarug,pbarvg,dtvr,p,pk,teta )
171c   ----------------------------------------------------------------
172c   Schema de Frederic Hourdin
173c   ----------------------------------------------------------------
174        else if(iadv(iq).eq.12) then
175c            Pas de temps adaptatif
176           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
177           if (n.GT.1) then
178           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
179     s             dtvr,'n=',n
180           endif
181           do indice=1,n
182            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
183           end do
184        else if(iadv(iq).eq.13) then
185c            Pas de temps adaptatif
186           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
187           if (n.GT.1) then
188           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
189     s             dtvr,'n=',n
190           endif
191          do indice=1,n
192            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
193          end do
194c   ----------------------------------------------------------------
195c   Schema de pente SLOPES
196c   ----------------------------------------------------------------
197        else if (iadv(iq).eq.20) then
198            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
199
200c   ----------------------------------------------------------------
201c   Schema de Prather
202c   ----------------------------------------------------------------
203        else if (iadv(iq).eq.30) then
204c            Pas de temps adaptatif
205           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
206           if (n.GT.1) then
207           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
208     s             dtvr,'n=',n
209           endif
210           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
211     s                     n,dtbon)
[960]212
[524]213c   ----------------------------------------------------------------
214c   Schemas PPM Lin et Rood
215c   ----------------------------------------------------------------
216         else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
217     s                     iadv(iq).LE.18)) then
218
219c        Test sur le flux horizontal
220c        Pas de temps adaptatif
221         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
222         if (n.GT.1) then
223         write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
224     s             dtvr,'n=',n
225         endif
226c        Test sur le flux vertical
227         CFLmaxz=0.
228         do l=2,llm
229           do ij=iip2,ip1jm
230            aaa=wg(ij,l)*dtvr/massem(ij,l)
231            CFLmaxz=max(CFLmaxz,aaa)
232            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
233            CFLmaxz=max(CFLmaxz,bbb)
234           enddo
235         enddo
236         if (CFLmaxz.GE.1) then
237            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
238         endif
239
240c-----------------------------------------------------------
241c        Ss-prg interface LMDZ.4->PPM3d
242c-----------------------------------------------------------
243
244          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
245     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
246     s                 unatppm,vnatppm,psppm)
247
248          do indice=1,n
249c---------------------------------------------------------------------
250c                         VL (version PPM) horiz. et PPM vert.
251c---------------------------------------------------------------------
252                if (iadv(iq).eq.11) then
253c                  Ss-prg PPM3d de Lin
254                  call ppm3d(1,qppm(1,1,iq),
255     s                       psppm,psppm,
256     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
257     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
258     s                       fill,dum,220.)
259
260c----------------------------------------------------------------------
261c                           Monotonic PPM
262c----------------------------------------------------------------------
263               else if (iadv(iq).eq.16) then
264c                  Ss-prg PPM3d de Lin
265                  call ppm3d(1,qppm(1,1,iq),
266     s                       psppm,psppm,
267     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
268     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
269     s                       fill,dum,220.)
270c---------------------------------------------------------------------
271
272c---------------------------------------------------------------------
273c                           Semi Monotonic PPM
274c---------------------------------------------------------------------
275               else if (iadv(iq).eq.17) then
276c                  Ss-prg PPM3d de Lin
277                  call ppm3d(1,qppm(1,1,iq),
278     s                       psppm,psppm,
279     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
280     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
281     s                       fill,dum,220.)
282c---------------------------------------------------------------------
283
284c---------------------------------------------------------------------
285c                         Positive Definite PPM
286c---------------------------------------------------------------------
287                else if (iadv(iq).eq.18) then
288c                  Ss-prg PPM3d de Lin
289                  call ppm3d(1,qppm(1,1,iq),
290     s                       psppm,psppm,
291     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
292     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
293     s                       fill,dum,220.)
294c---------------------------------------------------------------------
295                endif
296            enddo
297c-----------------------------------------------------------------
298c               Ss-prg interface PPM3d-LMDZ.4
299c-----------------------------------------------------------------
300                  call interpost(q(1,1,iq),qppm(1,1,iq))
301            endif
302c----------------------------------------------------------------------
303
304c-----------------------------------------------------------------
305c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
306c et Nord j=1
307c-----------------------------------------------------------------
308
309c                  call traceurpole(q(1,1,iq),massem)
310
311c calcul du temps cpu pour un schema donne
312
313c                  call clock(t_final)
[566]314cym                  tps_cpu=t_final-t_initial
315cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
[524]316
317       end DO
318
319
320c------------------------------------------------------------------
321c   on reinitialise a zero les flux de masse cumules
322c---------------------------------------------------
323          iadvtr=0
324
325       ENDIF ! if iadvtr.EQ.iapp_tracvl
326
327       RETURN
328       END
329
Note: See TracBrowser for help on using the repository browser.