source: LMDZ4/branches/pre_V3/libf/dyn3dpar/advtrac.F @ 789

Last change on this file since 789 was 630, checked in by Laurent Fairhead, 20 years ago

Import d'une version parallele de la dynamique YM
LF

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