source: LMDZ4/branches/LMDZ4V5.0-LF/libf/dyn3d/advtrac.F @ 4082

Last change on this file since 4082 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

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