source: LMDZ4/trunk/libf/dyn3d/advtrac.F @ 616

Last change on this file since 616 was 616, checked in by lmdzadmin, 19 years ago

Mise a jour pour INCA.2.0 Anne C
MAFi+LF

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