source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/advtrac.F @ 1183

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

Modif pour la compilation avec fcm sur Vargas SD
Corrections relatives au controle de histins et de l'utilisation des
options INST(X) pour l'appel au routines de IOIPSL FH
Mise en place de diagnostics sur les critères CFL pour l'advection
de traceurs. FH
Controle dans les .def de la dépendance verticale de l'efficacite
de la diffusion. Actif pour le moment uniquement avec ok_strato=y et llm=39. FH
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.7 KB
RevLine 
[524]1!
[1176]2! $Id: advtrac.F 1176 2009-06-11 08:54:10Z jghattas $
[524]3!
4c
5c
6      SUBROUTINE advtrac(pbaru,pbarv ,
7     *                   p,  masse,q,iapptrac,teta,
8     *                  flxw,
[762]9     *                  pk)
[524]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
[1114]17      USE infotrac
18
[524]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"
[1176]32#include "iniprint.h"
[524]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)
[1114]43      REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
[524]44      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
45      REAL pk(ip1jmp1,llm)
[960]46      REAL flxw(ip1jmp1,llm)
[524]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)
[1114]70      REAL qppm(iim*jjp1,llm,nqtot)
[524]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
[1176]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)
[524]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
[1176]147
[524]148c-------------------------------------------------------------------
[1176]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-------------------------------------------------------------------
[524]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-----------------------------------------------------------
[1114]228      do iq=1,nqtot
[524]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)
[960]287
[524]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)
[566]389cym                  tps_cpu=t_final-t_initial
390cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
[524]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.