source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3d/advtrac.F @ 5394

Last change on this file since 5394 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.7 KB
Line 
1!
2! $Id: advtrac.F 1279 2009-12-10 09:02:56Z evignon $
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
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#include "iniprint.h"
33
34c-------------------------------------------------------------------
35c     Arguments
36c-------------------------------------------------------------------
37c     Ajout PPM
38c--------------------------------------------------------
39      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
40c--------------------------------------------------------
41      INTEGER iapptrac
42      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
43      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
44      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
45      REAL pk(ip1jmp1,llm)
46      REAL flxw(ip1jmp1,llm)
47
48c-------------------------------------------------------------
49c     Variables locales
50c-------------------------------------------------------------
51
52      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
53      REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
54      REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
55      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
56      INTEGER iadvtr
57      INTEGER ij,l,iq,iiq
58      REAL zdpmin, zdpmax
59      EXTERNAL  minmax
60      SAVE iadvtr, massem, pbaruc, pbarvc
61      DATA iadvtr/0/
62c----------------------------------------------------------
63c     Rajouts pour PPM
64c----------------------------------------------------------
65      INTEGER indice,n
66      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
67      REAL CFLmaxz,aaa,bbb ! CFL maximum
68      REAL psppm(iim,jjp1) ! pression  au sol
69      REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
70      REAL qppm(iim*jjp1,llm,nqtot)
71      REAL fluxwppm(iim,jjp1,llm)
72      REAL apppm(llmp1), bpppm(llmp1)
73      LOGICAL dum,fill
74      DATA fill/.true./
75      DATA dum/.true./
76
77      integer,save :: countcfl=0
78      real cflx(ip1jmp1,llm)
79      real cfly(ip1jm,llm)
80      real cflz(ip1jmp1,llm)
81      real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
82
83      IF(iadvtr.EQ.0) THEN
84         CALL initial0(ijp1llm,pbaruc)
85         CALL initial0(ijmllm,pbarvc)
86      ENDIF
87
88c   accumulation des flux de masse horizontaux
89      DO l=1,llm
90         DO ij = 1,ip1jmp1
91            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
92         ENDDO
93         DO ij = 1,ip1jm
94            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
95         ENDDO
96      ENDDO
97
98c   selection de la masse instantannee des mailles avant le transport.
99      IF(iadvtr.EQ.0) THEN
100
101         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
102ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
103c
104      ENDIF
105
106      iadvtr   = iadvtr+1
107      iapptrac = iadvtr
108
109
110c   Test pour savoir si on advecte a ce pas de temps
111      IF ( iadvtr.EQ.iapp_tracvl ) THEN
112
113cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
114cc
115
116c   traitement des flux de masse avant advection.
117c     1. calcul de w
118c     2. groupement des mailles pres du pole.
119
120        CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
121
122      ! ... Flux de masse diaganostiques traceurs
123      flxw = wg / FLOAT(iapp_tracvl)
124
125c  test sur l'eventuelle creation de valeurs negatives de la masse
126         DO l=1,llm-1
127            DO ij = iip2+1,ip1jm
128              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
129     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
130     s                  +       wg(ij,l+1)  - wg(ij,l)
131            ENDDO
132            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
133            DO ij = iip2,ip1jm
134               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
135            ENDDO
136
137
138            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
139
140            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
141            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
142     s        '   MAX:', zdpmax
143            ENDIF
144
145         ENDDO
146
147
148c-------------------------------------------------------------------
149! Calcul des criteres CFL en X, Y et Z
150c-------------------------------------------------------------------
151
152      if (countcfl == 0. ) then
153          cflxmax(:)=0.
154          cflymax(:)=0.
155          cflzmax(:)=0.
156      endif
157
158      countcfl=countcfl+iapp_tracvl
159      cflx(:,:)=0.
160      cfly(:,:)=0.
161      cflz(:,:)=0.
162      do l=1,llm
163         do ij=iip2,ip1jm-1
164            if (pbarug(ij,l)>=0.) then
165                cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
166            else
167                cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
168            endif
169         enddo
170      enddo
171      do l=1,llm
172         do ij=iip2,ip1jm-1,iip1
173            cflx(ij+iip1,l)=cflx(ij,l)
174         enddo
175      enddo
176
177      do l=1,llm
178         do ij=1,ip1jm
179            if (pbarvg(ij,l)>=0.) then
180                cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
181            else
182                cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
183            endif
184         enddo
185      enddo
186
187      do l=2,llm
188         do ij=1,ip1jm
189            if (wg(ij,l)>=0.) then
190                cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
191            else
192                cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
193            endif
194         enddo
195      enddo
196
197      do l=1,llm
198         cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
199         cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
200         cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
201      enddo
202
203!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
204! Par defaut, on sort le diagnostic des CFL tous les jours.
205! Si on veut le sortir a chaque pas d'advection en cas de plantage
206!     if (countcfl==iapp_tracvl) then
207!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208      if (countcfl==day_step) then
209         do l=1,llm
210         write(lunout,*) 'L, CFLmax '
211     s   ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
212         enddo
213         countcfl=0
214      endif
215   
216c-------------------------------------------------------------------
217c   Advection proprement dite (Modification Le Croller (07/2001)
218c-------------------------------------------------------------------
219
220c----------------------------------------------------
221c        Calcul des moyennes basées sur la masse
222c----------------------------------------------------
223          call massbar(massem,massebx,masseby)         
224
225c-----------------------------------------------------------
226c     Appel des sous programmes d'advection
227c-----------------------------------------------------------
228      do iq=1,nqtot
229c        call clock(t_initial)
230        if(iadv(iq) == 0) cycle
231c   ----------------------------------------------------------------
232c   Schema de Van Leer I MUSCL
233c   ----------------------------------------------------------------
234        if(iadv(iq).eq.10) THEN
235            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
236
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.