source: LMDZ5/branches/LF-private/libf/dyn3d/advtrac.F90 @ 5361

Last change on this file since 5361 was 1817, checked in by lguez, 11 years ago

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