source: trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90 @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

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