source: LMDZ4/branches/V3_test/libf/dyn3d/advtrac.F @ 4082

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

Modifications pour rendre INCA plus independant de LMDZ ACo
LF

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