source: LMDZ6/trunk/libf/dyn3d/advtrac.F90 @ 4052

Last change on this file since 4052 was 4050, checked in by dcugnet, 3 years ago

Second commit for new tracers.

  • include most of the keys in the tracers descriptor vector "tracers(:)".
  • fix in phylmdiso/cv3_routines: fq_* variables were used where their fxt_* counterparts were expected.
  • multiple IF(nqdesc(iq)>0) and IF(nqfils(iq)>0) tests suppressed, because they are not needed: "do ... enddo" loops with 0 upper bound are not executed.
  • remove French accents from comments (encoding problem) in phylmdiso/cv3_routines and phylmdiso/cv30_routines.
  • modifications in "isotopes_verif_mod", where the call to function "iso_verif_tag17_q_deltad_chn" in "iso_verif_tag17_q_deltad_chn" was not detected at linking stage, although defined in the same module (?).
  • 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.4 KB
Line 
1! $Id: advtrac.F90 4050 2021-12-23 17:54:17Z dcugnet $
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, tracers, 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 "description.h"
22  include "iniprint.h"
23
24  !-------------------------------------------------------------------
25  !     Arguments
26  !-------------------------------------------------------------------
27  INTEGER,INTENT(OUT) :: iapptrac
28  REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
29  REAL,INTENT(IN) :: pbarv(ip1jm,llm)
30  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
31  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
32  REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
33  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
34  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
35  REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
36  !-------------------------------------------------------------------
37  !     Ajout PPM
38  !--------------------------------------------------------
39  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
40  !-------------------------------------------------------------
41  !     Variables locales
42  !-------------------------------------------------------------
43
44  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
45  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
46  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
47  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
48  INTEGER iadvtr
49  INTEGER ij,l,iq,iiq
50  REAL zdpmin, zdpmax
51  EXTERNAL  minmax
52  SAVE iadvtr, massem, pbaruc, pbarvc
53  DATA iadvtr/0/
54  !----------------------------------------------------------
55  !     Rajouts pour PPM
56  !----------------------------------------------------------
57  INTEGER indice,n
58  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
59  REAL CFLmaxz,aaa,bbb ! CFL maximum
60  REAL psppm(iim,jjp1) ! pression  au sol
61  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
62  REAL qppm(iim*jjp1,llm,nqtot)
63  REAL fluxwppm(iim,jjp1,llm)
64  REAL apppm(llmp1), bpppm(llmp1)
65  LOGICAL dum,fill
66  DATA fill/.true./
67  DATA dum/.true./
68
69  integer,save :: countcfl=0
70  real cflx(ip1jmp1,llm)
71  real cfly(ip1jm,llm)
72  real cflz(ip1jmp1,llm)
73  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
74  INTEGER :: iadv
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       iadv = tracers(iq)%iadv
230       SELECT CASE(iadv)
231         CASE(0); CYCLE
232         CASE(10)
233        !   ----------------------------------------------------------------
234        !   Schema de Van Leer I MUSCL
235        !   ----------------------------------------------------------------
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         CASE(14)
242           !   ----------------------------------------------------------------
243           !   Schema "pseudo amont" + test sur humidite specifique
244           !    pour la vapeur d'eau. F. Codron
245           !   ----------------------------------------------------------------
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         CASE(12)
252           !   ----------------------------------------------------------------
253           !   Schema de Frederic Hourdin
254           !   ----------------------------------------------------------------
255           !            Pas de temps adaptatif
256           call adaptdt(iadv,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         CASE(13)
265           !            Pas de temps adaptatif
266           call adaptdt(iadv,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         CASE(20)
275           !   ----------------------------------------------------------------
276           !   Schema de pente SLOPES
277           !   ----------------------------------------------------------------
278           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
279
280         CASE(30)
281           !   ----------------------------------------------------------------
282           !   Schema de Prather
283           !   ----------------------------------------------------------------
284           !            Pas de temps adaptatif
285           call adaptdt(iadv,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         CASE(11,16,17,18)
294           !   ----------------------------------------------------------------
295           !   Schemas PPM Lin et Rood
296           !   ----------------------------------------------------------------
297
298           !        Test sur le flux horizontal
299           !        Pas de temps adaptatif
300           call adaptdt(iadv,dtbon,n,pbarug,massem)
301           if (n.GT.1) then
302              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
303                   dtvr,'n=',n
304           endif
305           !        Test sur le flux vertical
306           CFLmaxz=0.
307           do l=2,llm
308              do ij=iip2,ip1jm
309                 aaa=wg(ij,l)*dtvr/massem(ij,l)
310                 CFLmaxz=max(CFLmaxz,aaa)
311                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
312                 CFLmaxz=max(CFLmaxz,bbb)
313              enddo
314           enddo
315           if (CFLmaxz.GE.1) then
316              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
317           endif
318
319           !-----------------------------------------------------------
320           !        Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
321           !-----------------------------------------------------------
322
323           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
324                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
325                unatppm,vnatppm,psppm)
326
327           do indice=1,n
328              !----------------------------------------------------------------
329              !                         VL (version PPM) horiz. et PPM vert.
330              !----------------------------------------------------------------
331             SELECT CASE(iadv)
332               CASE(11)
333                 call ppm3d(1,qppm(1,1,iq), &
334                      psppm,psppm, &
335                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
336                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
337                      fill,dum,220.)
338
339               CASE(16)
340                 !-------------------------------------------------------------
341                 !                           Monotonic PPM
342                 !-------------------------------------------------------------
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               CASE(17)
351                 !-------------------------------------------------------------
352                 !                           Semi Monotonic PPM
353                 !-------------------------------------------------------------
354                 call ppm3d(1,qppm(1,1,iq), &
355                      psppm,psppm, &
356                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
357                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
358                      fill,dum,220.)
359                 !-------------------------------------------------------------
360
361               CASE(18)
362                 !-------------------------------------------------------------
363                 !                         Positive Definite PPM
364                 !-------------------------------------------------------------
365                 call ppm3d(1,qppm(1,1,iq), &
366                      psppm,psppm, &
367                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
368                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
369                      fill,dum,220.)
370                 !-------------------------------------------------------------
371             END SELECT
372           enddo
373           !-----------------------------------------------------------------
374           !               Ss-prg interface PPM3d-LMDZ.4
375           !-----------------------------------------------------------------
376           call interpost(q(1,1,iq),qppm(1,1,iq))
377         END SELECT
378        !----------------------------------------------------------------------
379
380        !-----------------------------------------------------------------
381        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
382        ! et Nord j=1
383        !-----------------------------------------------------------------
384
385        !                  call traceurpole(q(1,1,iq),massem)
386
387        ! calcul du temps cpu pour un schema donne
388
389        !                  call clock(t_final)
390        !ym                  tps_cpu=t_final-t_initial
391        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
392
393     end DO
394
395     if (ok_iso_verif) then
396           write(*,*) 'advtrac 402'
397           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
398     endif !if (ok_iso_verif) then
399
400     !------------------------------------------------------------------
401     !   on reinitialise a zero les flux de masse cumules
402     !---------------------------------------------------
403     iadvtr=0
404
405  ENDIF ! if iadvtr.EQ.iapp_tracvl
406
407END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.