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

Last change on this file since 2600 was 2600, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.8 KB
Line 
1! $Id: advtrac.F90 2600 2016-07-23 05:45:38Z emillour $
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, ONLY: nqtot, iadv,nqperes,ok_iso_verif
12  USE control_mod, ONLY: iapp_tracvl, day_step
13  USE comconst_mod, ONLY: dtvr
14
15  IMPLICIT NONE
16  !
17  include "dimensions.h"
18  include "paramet.h"
19  include "comdissip.h"
20  include "comgeom2.h"
21  include "logic.h"
22  include "temps.h"
23  include "ener.h"
24  include "description.h"
25  include "iniprint.h"
26
27  !-------------------------------------------------------------------
28  !     Arguments
29  !-------------------------------------------------------------------
30  INTEGER,INTENT(OUT) :: iapptrac
31  REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
32  REAL,INTENT(IN) :: pbarv(ip1jm,llm)
33  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
34  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
35  REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
36  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
37  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
38  REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
39  !-------------------------------------------------------------------
40  !     Ajout PPM
41  !--------------------------------------------------------
42  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
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     pbaruc(:,:)=0
80     pbarvc(:,:)=0
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
224     if (ok_iso_verif) then
225           write(*,*) 'advtrac 227'
226           call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
227     endif !if (ok_iso_verif) then
228
229     do iq=1,nqperes
230        !        call clock(t_initial)
231        if(iadv(iq) == 0) cycle
232        !   ----------------------------------------------------------------
233        !   Schema de Van Leer I MUSCL
234        !   ----------------------------------------------------------------
235        if(iadv(iq).eq.10) THEN
236           ! CRisi: on fait passer tout q pour avoir acces aux fils
237           
238           !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
239           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
240           
241           !   ----------------------------------------------------------------
242           !   Schema "pseudo amont" + test sur humidite specifique
243           !    pour la vapeur d'eau. F. Codron
244           !   ----------------------------------------------------------------
245        else if(iadv(iq).eq.14) then
246           !
247           !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
248           CALL vlspltqs( q, 2., massem, wg , &
249                pbarug,pbarvg,dtvr,p,pk,teta,iq)
250           
251           !   ----------------------------------------------------------------
252           !   Schema de Frederic Hourdin
253           !   ----------------------------------------------------------------
254        else if(iadv(iq).eq.12) then
255           !            Pas de temps adaptatif
256           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
257           if (n.GT.1) then
258              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
259                   dtvr,'n=',n
260           endif
261           do indice=1,n
262              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
263           end do
264        else if(iadv(iq).eq.13) then
265           !            Pas de temps adaptatif
266           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
267           if (n.GT.1) then
268              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
269                   dtvr,'n=',n
270           endif
271           do indice=1,n
272              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
273           end do
274           !   ----------------------------------------------------------------
275           !   Schema de pente SLOPES
276           !   ----------------------------------------------------------------
277        else if (iadv(iq).eq.20) then
278           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
279
280           !   ----------------------------------------------------------------
281           !   Schema de Prather
282           !   ----------------------------------------------------------------
283        else if (iadv(iq).eq.30) then
284           !            Pas de temps adaptatif
285           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
286           if (n.GT.1) then
287              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
288                   dtvr,'n=',n
289           endif
290           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
291                n,dtbon)
292
293           !   ----------------------------------------------------------------
294           !   Schemas PPM Lin et Rood
295           !   ----------------------------------------------------------------
296        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
297             iadv(iq).LE.18)) then
298
299           !        Test sur le flux horizontal
300           !        Pas de temps adaptatif
301           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
302           if (n.GT.1) then
303              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
304                   dtvr,'n=',n
305           endif
306           !        Test sur le flux vertical
307           CFLmaxz=0.
308           do l=2,llm
309              do ij=iip2,ip1jm
310                 aaa=wg(ij,l)*dtvr/massem(ij,l)
311                 CFLmaxz=max(CFLmaxz,aaa)
312                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
313                 CFLmaxz=max(CFLmaxz,bbb)
314              enddo
315           enddo
316           if (CFLmaxz.GE.1) then
317              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
318           endif
319
320           !-----------------------------------------------------------
321           !        Ss-prg interface LMDZ.4->PPM3d
322           !-----------------------------------------------------------
323
324           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
325                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
326                unatppm,vnatppm,psppm)
327
328           do indice=1,n
329              !----------------------------------------------------------------
330              !                         VL (version PPM) horiz. et PPM vert.
331              !----------------------------------------------------------------
332              if (iadv(iq).eq.11) then
333                 !                  Ss-prg PPM3d de Lin
334                 call ppm3d(1,qppm(1,1,iq), &
335                      psppm,psppm, &
336                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
337                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
338                      fill,dum,220.)
339
340                 !-------------------------------------------------------------
341                 !                           Monotonic PPM
342                 !-------------------------------------------------------------
343              else if (iadv(iq).eq.16) then
344                 !                  Ss-prg PPM3d de Lin
345                 call ppm3d(1,qppm(1,1,iq), &
346                      psppm,psppm, &
347                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
348                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
349                      fill,dum,220.)
350                 !-------------------------------------------------------------
351
352                 !-------------------------------------------------------------
353                 !                           Semi Monotonic PPM
354                 !-------------------------------------------------------------
355              else if (iadv(iq).eq.17) then
356                 !                  Ss-prg PPM3d de Lin
357                 call ppm3d(1,qppm(1,1,iq), &
358                      psppm,psppm, &
359                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
360                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
361                      fill,dum,220.)
362                 !-------------------------------------------------------------
363
364                 !-------------------------------------------------------------
365                 !                         Positive Definite PPM
366                 !-------------------------------------------------------------
367              else if (iadv(iq).eq.18) then
368                 !                  Ss-prg PPM3d de Lin
369                 call ppm3d(1,qppm(1,1,iq), &
370                      psppm,psppm, &
371                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
372                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
373                      fill,dum,220.)
374                 !-------------------------------------------------------------
375              endif
376           enddo
377           !-----------------------------------------------------------------
378           !               Ss-prg interface PPM3d-LMDZ.4
379           !-----------------------------------------------------------------
380           call interpost(q(1,1,iq),qppm(1,1,iq))
381        endif
382        !----------------------------------------------------------------------
383
384        !-----------------------------------------------------------------
385        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
386        ! et Nord j=1
387        !-----------------------------------------------------------------
388
389        !                  call traceurpole(q(1,1,iq),massem)
390
391        ! calcul du temps cpu pour un schema donne
392
393        !                  call clock(t_final)
394        !ym                  tps_cpu=t_final-t_initial
395        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
396
397     end DO
398
399     if (ok_iso_verif) then
400           write(*,*) 'advtrac 402'
401           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
402     endif !if (ok_iso_verif) then
403
404     !------------------------------------------------------------------
405     !   on reinitialise a zero les flux de masse cumules
406     !---------------------------------------------------
407     iadvtr=0
408
409  ENDIF ! if iadvtr.EQ.iapp_tracvl
410
411END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.