source: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/advtrac.F @ 1818

Last change on this file since 1818 was 1454, checked in by Laurent Fairhead, 14 years ago

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434


Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

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