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

Last change on this file since 2436 was 1508, checked in by emillour, 10 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

  • in dynlonlat_phylonlat:
  • add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)
  • in misc:
  • follow up updates on wxios.F (add missing_val module variable)
  • in dyn3d_common:
  • pression.F => pression.F90
  • misc_mod.F90: moved from misc to dyn3d_common
  • added new iso_verif_dyn.F
  • covcont.F => covcont.F90
  • infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
  • dynetat0.F => dynetat0.F90 with some code factorization
  • dynredem.F => dynredem.F90 with some code factorization
  • added dynredem_mod.F90: routines used by dynredem
  • iniacademic.F90 : added isotopes-related initialization for Earth case
  • in dyn3d:
  • added check_isotopes.F
  • modified (isotopes) advtrac.F90, caladvtrac.F
  • guide_mod.F90: ported updates
  • leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
  • qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
  • vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
  • vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)
  • in dyn3dpar:
  • leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
  • dynredem_p.F => dynredem_p.F90 and some code factorization
  • and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

EM
File size: 14.7 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, nqperes, ok_iso_verif
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     pbaruc(:,:)=0
77     pbarvc(:,:)=0
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     if (ok_iso_verif) then
221           write(*,*) 'advtrac 227'
222           call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
223     endif !if (ok_iso_verif) then
224
225     do iq=1,nqperes
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           ! CRisi: on fait passer tout q pour avoir acces aux fils
233           
234           !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
235           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
236
237           !   ----------------------------------------------------------------
238           !   Schema "pseudo amont" + test sur humidite specifique
239           !    pour la vapeur d'eau. F. Codron
240           !   ----------------------------------------------------------------
241        else if(iadv(iq).eq.14) then
242           !
243           !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
244           CALL vlspltqs( q, 2., massem, wg , &
245                pbarug,pbarvg,dtvr,p,pk,teta,iq)
246
247           !   ----------------------------------------------------------------
248           !   Schema de Frederic Hourdin
249           !   ----------------------------------------------------------------
250        else if(iadv(iq).eq.12) 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,1)
259           end do
260        else if(iadv(iq).eq.13) then
261           !            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                   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
270           !   ----------------------------------------------------------------
271           !   Schema de pente SLOPES
272           !   ----------------------------------------------------------------
273        else if (iadv(iq).eq.20) then
274           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
275
276           !   ----------------------------------------------------------------
277           !   Schema de Prather
278           !   ----------------------------------------------------------------
279        else if (iadv(iq).eq.30) then
280           !            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                   dtvr,'n=',n
285           endif
286           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
287                n,dtbon)
288
289           !   ----------------------------------------------------------------
290           !   Schemas PPM Lin et Rood
291           !   ----------------------------------------------------------------
292        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
293             iadv(iq).LE.18)) then
294
295           !        Test sur le flux horizontal
296           !        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                   dtvr,'n=',n
301           endif
302           !        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
316           !-----------------------------------------------------------
317           !        Ss-prg interface LMDZ.4->PPM3d
318           !-----------------------------------------------------------
319
320           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
321                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
322                unatppm,vnatppm,psppm)
323
324           do indice=1,n
325              !----------------------------------------------------------------
326              !                         VL (version PPM) horiz. et PPM vert.
327              !----------------------------------------------------------------
328              if (iadv(iq).eq.11) then
329                 !                  Ss-prg PPM3d de Lin
330                 call ppm3d(1,qppm(1,1,iq), &
331                      psppm,psppm, &
332                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
333                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
334                      fill,dum,220.)
335
336                 !-------------------------------------------------------------
337                 !                           Monotonic PPM
338                 !-------------------------------------------------------------
339              else if (iadv(iq).eq.16) then
340                 !                  Ss-prg PPM3d de Lin
341                 call ppm3d(1,qppm(1,1,iq), &
342                      psppm,psppm, &
343                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
344                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
345                      fill,dum,220.)
346                 !-------------------------------------------------------------
347
348                 !-------------------------------------------------------------
349                 !                           Semi Monotonic PPM
350                 !-------------------------------------------------------------
351              else if (iadv(iq).eq.17) then
352                 !                  Ss-prg PPM3d de Lin
353                 call ppm3d(1,qppm(1,1,iq), &
354                      psppm,psppm, &
355                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
356                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
357                      fill,dum,220.)
358                 !-------------------------------------------------------------
359
360                 !-------------------------------------------------------------
361                 !                         Positive Definite PPM
362                 !-------------------------------------------------------------
363              else if (iadv(iq).eq.18) then
364                 !                  Ss-prg PPM3d de Lin
365                 call ppm3d(1,qppm(1,1,iq), &
366                      psppm,psppm, &
367                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
368                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
369                      fill,dum,220.)
370                 !-------------------------------------------------------------
371              endif
372           enddo
373           !-----------------------------------------------------------------
374           !               Ss-prg interface PPM3d-LMDZ.4
375           !-----------------------------------------------------------------
376           call interpost(q(1,1,iq),qppm(1,1,iq))
377        endif
378        !----------------------------------------------------------------------
379
380        !-----------------------------------------------------------------
381        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
382        ! et Nord j=1
383        !-----------------------------------------------------------------
384
385        !                  call traceurpole(q(1,1,iq),massem)
386
387        ! calcul du temps cpu pour un schema donne
388
389        !                  call clock(t_final)
390        !ym                  tps_cpu=t_final-t_initial
391        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
392
393     end DO
394
395     if (ok_iso_verif) then
396           write(*,*) 'advtrac 402'
397           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
398     endif !if (ok_iso_verif) then
399
400     !------------------------------------------------------------------
401     !   on reinitialise a zero les flux de masse cumules
402     !---------------------------------------------------
403     iadvtr=0
404
405  ENDIF ! if iadvtr.EQ.iapp_tracvl
406
407END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.