source: LMDZ5/trunk/libf/dyn3d/advtrac.F90 @ 1839

Last change on this file since 1839 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
RevLine 
[1279]1! $Id: advtrac.F90 1817 2013-07-25 12:49:07Z fairhead $
[1146]2
[1549]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
[524]13
14
[1549]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"
[524]28
[1549]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)
[524]42
[1549]43  !-------------------------------------------------------------
44  !     Variables locales
45  !-------------------------------------------------------------
[524]46
[1549]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./
[524]71
[1549]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)
[524]77
[1549]78  IF(iadvtr.EQ.0) THEN
79     CALL initial0(ijp1llm,pbaruc)
80     CALL initial0(ijmllm,pbarvc)
81  ENDIF
[524]82
[1549]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
[524]92
[1549]93  !   selection de la masse instantannee des mailles avant le transport.
94  IF(iadvtr.EQ.0) THEN
[524]95
[1549]96     CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
97     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
98     !
99  ENDIF
[524]100
[1549]101  iadvtr   = iadvtr+1
102  iapptrac = iadvtr
[524]103
104
[1549]105  !   Test pour savoir si on advecte a ce pas de temps
106  IF ( iadvtr.EQ.iapp_tracvl ) THEN
[524]107
[1549]108     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
109     !c
[524]110
[1549]111     !   traitement des flux de masse avant advection.
112     !     1. calcul de w
113     !     2. groupement des mailles pres du pole.
[524]114
[1549]115     CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
[524]116
[1549]117     ! ... Flux de masse diaganostiques traceurs
118     flxw = wg / REAL(iapp_tracvl)
[524]119
[1549]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
[524]131
132
[1549]133        CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
[524]134
[1549]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
[1279]139
[1549]140     ENDDO
[1279]141
142
[1549]143     !-------------------------------------------------------------------
144     ! Calcul des criteres CFL en X, Y et Z
145     !-------------------------------------------------------------------
[1279]146
[1549]147     if (countcfl == 0. ) then
148        cflxmax(:)=0.
149        cflymax(:)=0.
150        cflzmax(:)=0.
151     endif
[1279]152
[1549]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
[1279]171
[1549]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
[1279]181
[1549]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
[1279]198!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1549]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
[1279]202!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1549]203     if (countcfl==day_step) then
204        do l=1,llm
[1817]205           write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
206                cflzmax(l)
[1549]207        enddo
208        countcfl=0
209     endif
[524]210
[1549]211     !-------------------------------------------------------------------
212     !   Advection proprement dite (Modification Le Croller (07/2001)
213     !-------------------------------------------------------------------
[524]214
[1549]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)
[524]225        if(iadv(iq) == 0) cycle
[1549]226        !   ----------------------------------------------------------------
227        !   Schema de Van Leer I MUSCL
228        !   ----------------------------------------------------------------
[524]229        if(iadv(iq).eq.10) THEN
[1549]230           call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[524]231
[1549]232           !   ----------------------------------------------------------------
233           !   Schema "pseudo amont" + test sur humidite specifique
234           !    pour la vapeur d'eau. F. Codron
235           !   ----------------------------------------------------------------
[524]236        else if(iadv(iq).eq.14) then
[1549]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           !   ----------------------------------------------------------------
[524]243        else if(iadv(iq).eq.12) then
[1549]244           !            Pas de temps adaptatif
[524]245           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
246           if (n.GT.1) then
[1549]247              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
248                   dtvr,'n=',n
[524]249           endif
250           do indice=1,n
[1549]251              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
[524]252           end do
253        else if(iadv(iq).eq.13) then
[1549]254           !            Pas de temps adaptatif
[524]255           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
256           if (n.GT.1) then
[1549]257              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
258                   dtvr,'n=',n
[524]259           endif
[1549]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           !   ----------------------------------------------------------------
[524]266        else if (iadv(iq).eq.20) then
[1549]267           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[524]268
[1549]269           !   ----------------------------------------------------------------
270           !   Schema de Prather
271           !   ----------------------------------------------------------------
[524]272        else if (iadv(iq).eq.30) then
[1549]273           !            Pas de temps adaptatif
[524]274           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
275           if (n.GT.1) then
[1549]276              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
277                   dtvr,'n=',n
[524]278           endif
[1549]279           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
280                n,dtbon)
[960]281
[1549]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
[524]287
[1549]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
[524]304           enddo
[1549]305           if (CFLmaxz.GE.1) then
306              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
307           endif
[524]308
[1549]309           !-----------------------------------------------------------
310           !        Ss-prg interface LMDZ.4->PPM3d
311           !-----------------------------------------------------------
[524]312
[1549]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)
[524]316
[1549]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.)
[524]328
[1549]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                 !-------------------------------------------------------------
[524]340
[1549]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                 !-------------------------------------------------------------
[524]352
[1549]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        !----------------------------------------------------------------------
[524]372
[1549]373        !-----------------------------------------------------------------
374        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
375        ! et Nord j=1
376        !-----------------------------------------------------------------
[524]377
[1549]378        !                  call traceurpole(q(1,1,iq),massem)
[524]379
[1549]380        ! calcul du temps cpu pour un schema donne
[524]381
[1549]382        !                  call clock(t_final)
383        !ym                  tps_cpu=t_final-t_initial
384        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
[524]385
[1549]386     end DO
[524]387
388
[1549]389     !------------------------------------------------------------------
390     !   on reinitialise a zero les flux de masse cumules
391     !---------------------------------------------------
392     iadvtr=0
[524]393
[1549]394  ENDIF ! if iadvtr.EQ.iapp_tracvl
[524]395
[1549]396END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.