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

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

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