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

Last change on this file since 2146 was 1987, checked in by Ehouarn Millour, 11 years ago

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
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.3 KB
Line 
1! $Id: advtrac.F90 1987 2014-02-24 15:05:47Z idelkadi $
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
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     CALL initial0(ijp1llm,pbaruc)
82     CALL initial0(ijmllm,pbarvc)
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     do iq=1,nqtot
226        !        call clock(t_initial)
227        if(iadv(iq) == 0) cycle
228        !   ----------------------------------------------------------------
229        !   Schema de Van Leer I MUSCL
230        !   ----------------------------------------------------------------
231        if(iadv(iq).eq.10) THEN
232           call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
233
234           !   ----------------------------------------------------------------
235           !   Schema "pseudo amont" + test sur humidite specifique
236           !    pour la vapeur d'eau. F. Codron
237           !   ----------------------------------------------------------------
238        else if(iadv(iq).eq.14) then
239           !
240           CALL vlspltqs( q(1,1,1), 2., massem, wg , &
241                pbarug,pbarvg,dtvr,p,pk,teta )
242           !   ----------------------------------------------------------------
243           !   Schema de Frederic Hourdin
244           !   ----------------------------------------------------------------
245        else if(iadv(iq).eq.12) then
246           !            Pas de temps adaptatif
247           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
248           if (n.GT.1) then
249              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
250                   dtvr,'n=',n
251           endif
252           do indice=1,n
253              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
254           end do
255        else if(iadv(iq).eq.13) 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,2)
264           end do
265           !   ----------------------------------------------------------------
266           !   Schema de pente SLOPES
267           !   ----------------------------------------------------------------
268        else if (iadv(iq).eq.20) then
269           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
270
271           !   ----------------------------------------------------------------
272           !   Schema de Prather
273           !   ----------------------------------------------------------------
274        else if (iadv(iq).eq.30) then
275           !            Pas de temps adaptatif
276           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
277           if (n.GT.1) then
278              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
279                   dtvr,'n=',n
280           endif
281           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
282                n,dtbon)
283
284           !   ----------------------------------------------------------------
285           !   Schemas PPM Lin et Rood
286           !   ----------------------------------------------------------------
287        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
288             iadv(iq).LE.18)) then
289
290           !        Test sur le flux horizontal
291           !        Pas de temps adaptatif
292           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
293           if (n.GT.1) then
294              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
295                   dtvr,'n=',n
296           endif
297           !        Test sur le flux vertical
298           CFLmaxz=0.
299           do l=2,llm
300              do ij=iip2,ip1jm
301                 aaa=wg(ij,l)*dtvr/massem(ij,l)
302                 CFLmaxz=max(CFLmaxz,aaa)
303                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
304                 CFLmaxz=max(CFLmaxz,bbb)
305              enddo
306           enddo
307           if (CFLmaxz.GE.1) then
308              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
309           endif
310
311           !-----------------------------------------------------------
312           !        Ss-prg interface LMDZ.4->PPM3d
313           !-----------------------------------------------------------
314
315           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
316                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
317                unatppm,vnatppm,psppm)
318
319           do indice=1,n
320              !----------------------------------------------------------------
321              !                         VL (version PPM) horiz. et PPM vert.
322              !----------------------------------------------------------------
323              if (iadv(iq).eq.11) then
324                 !                  Ss-prg PPM3d de Lin
325                 call ppm3d(1,qppm(1,1,iq), &
326                      psppm,psppm, &
327                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
328                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
329                      fill,dum,220.)
330
331                 !-------------------------------------------------------------
332                 !                           Monotonic PPM
333                 !-------------------------------------------------------------
334              else if (iadv(iq).eq.16) then
335                 !                  Ss-prg PPM3d de Lin
336                 call ppm3d(1,qppm(1,1,iq), &
337                      psppm,psppm, &
338                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
339                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
340                      fill,dum,220.)
341                 !-------------------------------------------------------------
342
343                 !-------------------------------------------------------------
344                 !                           Semi Monotonic PPM
345                 !-------------------------------------------------------------
346              else if (iadv(iq).eq.17) then
347                 !                  Ss-prg PPM3d de Lin
348                 call ppm3d(1,qppm(1,1,iq), &
349                      psppm,psppm, &
350                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
351                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
352                      fill,dum,220.)
353                 !-------------------------------------------------------------
354
355                 !-------------------------------------------------------------
356                 !                         Positive Definite PPM
357                 !-------------------------------------------------------------
358              else if (iadv(iq).eq.18) then
359                 !                  Ss-prg PPM3d de Lin
360                 call ppm3d(1,qppm(1,1,iq), &
361                      psppm,psppm, &
362                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
363                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
364                      fill,dum,220.)
365                 !-------------------------------------------------------------
366              endif
367           enddo
368           !-----------------------------------------------------------------
369           !               Ss-prg interface PPM3d-LMDZ.4
370           !-----------------------------------------------------------------
371           call interpost(q(1,1,iq),qppm(1,1,iq))
372        endif
373        !----------------------------------------------------------------------
374
375        !-----------------------------------------------------------------
376        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
377        ! et Nord j=1
378        !-----------------------------------------------------------------
379
380        !                  call traceurpole(q(1,1,iq),massem)
381
382        ! calcul du temps cpu pour un schema donne
383
384        !                  call clock(t_final)
385        !ym                  tps_cpu=t_final-t_initial
386        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
387
388     end DO
389
390
391     !------------------------------------------------------------------
392     !   on reinitialise a zero les flux de masse cumules
393     !---------------------------------------------------
394     iadvtr=0
395
396  ENDIF ! if iadvtr.EQ.iapp_tracvl
397
398END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.