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

Last change on this file since 2279 was 2270, checked in by crisi, 10 years ago

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi

  • 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.7 KB
Line 
1! $Id: advtrac.F90 2270 2015-05-07 15:45:04Z musat $
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           call check_isotopes_seq(q,1,ip1jmp1,'advtrac 162')
228     endif !if (ok_iso_verif) then
229
230     do iq=1,nqperes
231        !        call clock(t_initial)
232        if(iadv(iq) == 0) cycle
233        !   ----------------------------------------------------------------
234        !   Schema de Van Leer I MUSCL
235        !   ----------------------------------------------------------------
236        if(iadv(iq).eq.10) THEN
237           ! CRisi: on fait passer tout q pour avoir acces aux fils
238           
239           write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
240           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
241           
242           !   ----------------------------------------------------------------
243           !   Schema "pseudo amont" + test sur humidite specifique
244           !    pour la vapeur d'eau. F. Codron
245           !   ----------------------------------------------------------------
246        else if(iadv(iq).eq.14) then
247           !
248           write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
249           CALL vlspltqs( q, 2., massem, wg , &
250                pbarug,pbarvg,dtvr,p,pk,teta,iq)
251           
252           !   ----------------------------------------------------------------
253           !   Schema de Frederic Hourdin
254           !   ----------------------------------------------------------------
255        else if(iadv(iq).eq.12) then
256           !            Pas de temps adaptatif
257           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
258           if (n.GT.1) then
259              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
260                   dtvr,'n=',n
261           endif
262           do indice=1,n
263              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
264           end do
265        else if(iadv(iq).eq.13) then
266           !            Pas de temps adaptatif
267           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
268           if (n.GT.1) then
269              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
270                   dtvr,'n=',n
271           endif
272           do indice=1,n
273              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
274           end do
275           !   ----------------------------------------------------------------
276           !   Schema de pente SLOPES
277           !   ----------------------------------------------------------------
278        else if (iadv(iq).eq.20) then
279           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
280
281           !   ----------------------------------------------------------------
282           !   Schema de Prather
283           !   ----------------------------------------------------------------
284        else if (iadv(iq).eq.30) then
285           !            Pas de temps adaptatif
286           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
287           if (n.GT.1) then
288              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
289                   dtvr,'n=',n
290           endif
291           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
292                n,dtbon)
293
294           !   ----------------------------------------------------------------
295           !   Schemas PPM Lin et Rood
296           !   ----------------------------------------------------------------
297        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
298             iadv(iq).LE.18)) then
299
300           !        Test sur le flux horizontal
301           !        Pas de temps adaptatif
302           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
303           if (n.GT.1) then
304              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
305                   dtvr,'n=',n
306           endif
307           !        Test sur le flux vertical
308           CFLmaxz=0.
309           do l=2,llm
310              do ij=iip2,ip1jm
311                 aaa=wg(ij,l)*dtvr/massem(ij,l)
312                 CFLmaxz=max(CFLmaxz,aaa)
313                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
314                 CFLmaxz=max(CFLmaxz,bbb)
315              enddo
316           enddo
317           if (CFLmaxz.GE.1) then
318              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
319           endif
320
321           !-----------------------------------------------------------
322           !        Ss-prg interface LMDZ.4->PPM3d
323           !-----------------------------------------------------------
324
325           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
326                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
327                unatppm,vnatppm,psppm)
328
329           do indice=1,n
330              !----------------------------------------------------------------
331              !                         VL (version PPM) horiz. et PPM vert.
332              !----------------------------------------------------------------
333              if (iadv(iq).eq.11) then
334                 !                  Ss-prg PPM3d de Lin
335                 call ppm3d(1,qppm(1,1,iq), &
336                      psppm,psppm, &
337                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
338                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
339                      fill,dum,220.)
340
341                 !-------------------------------------------------------------
342                 !                           Monotonic PPM
343                 !-------------------------------------------------------------
344              else if (iadv(iq).eq.16) then
345                 !                  Ss-prg PPM3d de Lin
346                 call ppm3d(1,qppm(1,1,iq), &
347                      psppm,psppm, &
348                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
349                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
350                      fill,dum,220.)
351                 !-------------------------------------------------------------
352
353                 !-------------------------------------------------------------
354                 !                           Semi Monotonic PPM
355                 !-------------------------------------------------------------
356              else if (iadv(iq).eq.17) then
357                 !                  Ss-prg PPM3d de Lin
358                 call ppm3d(1,qppm(1,1,iq), &
359                      psppm,psppm, &
360                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
361                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
362                      fill,dum,220.)
363                 !-------------------------------------------------------------
364
365                 !-------------------------------------------------------------
366                 !                         Positive Definite PPM
367                 !-------------------------------------------------------------
368              else if (iadv(iq).eq.18) then
369                 !                  Ss-prg PPM3d de Lin
370                 call ppm3d(1,qppm(1,1,iq), &
371                      psppm,psppm, &
372                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
373                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
374                      fill,dum,220.)
375                 !-------------------------------------------------------------
376              endif
377           enddo
378           !-----------------------------------------------------------------
379           !               Ss-prg interface PPM3d-LMDZ.4
380           !-----------------------------------------------------------------
381           call interpost(q(1,1,iq),qppm(1,1,iq))
382        endif
383        !----------------------------------------------------------------------
384
385        !-----------------------------------------------------------------
386        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
387        ! et Nord j=1
388        !-----------------------------------------------------------------
389
390        !                  call traceurpole(q(1,1,iq),massem)
391
392        ! calcul du temps cpu pour un schema donne
393
394        !                  call clock(t_final)
395        !ym                  tps_cpu=t_final-t_initial
396        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
397
398     end DO
399
400     if (ok_iso_verif) then
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.