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

Last change on this file since 2286 was 2286, checked in by Ehouarn Millour, 9 years ago
  • fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
  • comment out many invasive debug writes.

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