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

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

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