LMDZ
advtrac.F90
Go to the documentation of this file.
1 ! $Id: advtrac.F90 2286 2015-05-20 13:27:07Z emillour $
2 
3 SUBROUTINE 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 
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  pbaruc(:,:)=0
82  pbarvc(:,:)=0
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 
226  if (ok_iso_verif) then
227  write(*,*) 'advtrac 227'
228  call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
229  endif !if (ok_iso_verif) then
230 
231  do iq=1,nqperes
232  ! call clock(t_initial)
233  if(iadv(iq) == 0) cycle
234  ! ----------------------------------------------------------------
235  ! Schema de Van Leer I MUSCL
236  ! ----------------------------------------------------------------
237  if(iadv(iq).eq.10) THEN
238  ! CRisi: on fait passer tout q pour avoir acces aux fils
239 
240  !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
241  call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
242 
243  ! ----------------------------------------------------------------
244  ! Schema "pseudo amont" + test sur humidite specifique
245  ! pour la vapeur d'eau. F. Codron
246  ! ----------------------------------------------------------------
247  else if(iadv(iq).eq.14) then
248  !
249  !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
250  CALL vlspltqs( q, 2., massem, wg , &
251  pbarug,pbarvg,dtvr,p,pk,teta,iq)
252 
253  ! ----------------------------------------------------------------
254  ! Schema de Frederic Hourdin
255  ! ----------------------------------------------------------------
256  else if(iadv(iq).eq.12) then
257  ! Pas de temps adaptatif
258  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
259  if (n.GT.1) then
260  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
261  dtvr,'n=',n
262  endif
263  do indice=1,n
264  call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
265  end do
266  else if(iadv(iq).eq.13) then
267  ! Pas de temps adaptatif
268  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
269  if (n.GT.1) then
270  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
271  dtvr,'n=',n
272  endif
273  do indice=1,n
274  call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
275  end do
276  ! ----------------------------------------------------------------
277  ! Schema de pente SLOPES
278  ! ----------------------------------------------------------------
279  else if (iadv(iq).eq.20) then
280  call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
281 
282  ! ----------------------------------------------------------------
283  ! Schema de Prather
284  ! ----------------------------------------------------------------
285  else if (iadv(iq).eq.30) then
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  call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
293  n,dtbon)
294 
295  ! ----------------------------------------------------------------
296  ! Schemas PPM Lin et Rood
297  ! ----------------------------------------------------------------
298  else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
299  iadv(iq).LE.18)) then
300 
301  ! Test sur le flux horizontal
302  ! Pas de temps adaptatif
303  call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
304  if (n.GT.1) then
305  write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
306  dtvr,'n=',n
307  endif
308  ! Test sur le flux vertical
309  cflmaxz=0.
310  do l=2,llm
311  do ij=iip2,ip1jm
312  aaa=wg(ij,l)*dtvr/massem(ij,l)
313  cflmaxz=max(cflmaxz,aaa)
314  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
315  cflmaxz=max(cflmaxz,bbb)
316  enddo
317  enddo
318  if (cflmaxz.GE.1) then
319  write(*,*) 'WARNING vertical','CFLmaxz=', cflmaxz
320  endif
321 
322  !-----------------------------------------------------------
323  ! Ss-prg interface LMDZ.4->PPM3d
324  !-----------------------------------------------------------
325 
326  call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
327  apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
328  unatppm,vnatppm,psppm)
329 
330  do indice=1,n
331  !----------------------------------------------------------------
332  ! VL (version PPM) horiz. et PPM vert.
333  !----------------------------------------------------------------
334  if (iadv(iq).eq.11) then
335  ! Ss-prg PPM3d de Lin
336  call ppm3d(1,qppm(1,1,iq), &
337  psppm,psppm, &
338  unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
339  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
340  fill,dum,220.)
341 
342  !-------------------------------------------------------------
343  ! Monotonic PPM
344  !-------------------------------------------------------------
345  else if (iadv(iq).eq.16) then
346  ! Ss-prg PPM3d de Lin
347  call ppm3d(1,qppm(1,1,iq), &
348  psppm,psppm, &
349  unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
350  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
351  fill,dum,220.)
352  !-------------------------------------------------------------
353 
354  !-------------------------------------------------------------
355  ! Semi Monotonic PPM
356  !-------------------------------------------------------------
357  else if (iadv(iq).eq.17) then
358  ! Ss-prg PPM3d de Lin
359  call ppm3d(1,qppm(1,1,iq), &
360  psppm,psppm, &
361  unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
362  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
363  fill,dum,220.)
364  !-------------------------------------------------------------
365 
366  !-------------------------------------------------------------
367  ! Positive Definite PPM
368  !-------------------------------------------------------------
369  else if (iadv(iq).eq.18) then
370  ! Ss-prg PPM3d de Lin
371  call ppm3d(1,qppm(1,1,iq), &
372  psppm,psppm, &
373  unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
374  iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
375  fill,dum,220.)
376  !-------------------------------------------------------------
377  endif
378  enddo
379  !-----------------------------------------------------------------
380  ! Ss-prg interface PPM3d-LMDZ.4
381  !-----------------------------------------------------------------
382  call interpost(q(1,1,iq),qppm(1,1,iq))
383  endif
384  !----------------------------------------------------------------------
385 
386  !-----------------------------------------------------------------
387  ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
388  ! et Nord j=1
389  !-----------------------------------------------------------------
390 
391  ! call traceurpole(q(1,1,iq),massem)
392 
393  ! calcul du temps cpu pour un schema donne
394 
395  ! call clock(t_final)
396  !ym tps_cpu=t_final-t_initial
397  !ym cpuadv(iq)=cpuadv(iq)+tps_cpu
398 
399  end DO
400 
401  if (ok_iso_verif) then
402  write(*,*) 'advtrac 402'
403  call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
404  endif !if (ok_iso_verif) then
405 
406  !------------------------------------------------------------------
407  ! on reinitialise a zero les flux de masse cumules
408  !---------------------------------------------------
409  iadvtr=0
410 
411  ENDIF ! if iadvtr.EQ.iapp_tracvl
412 
413 END SUBROUTINE advtrac
pentes_ini
subroutine pentes_ini(q, w, masse, pbaru, pbarv, mode)
Definition: pentes_ini.F:5
control_mod::iapp_tracvl
integer, save iapp_tracvl
Definition: control_mod.F90:17
iip2
!$Header iip2
Definition: paramet.h:14
ip1jmp1
!$Header llmm1 INTEGER ip1jmp1
Definition: paramet.h:14
interpost
subroutine interpost(q, qppm)
Definition: interpost.F:5
adaptdt
subroutine adaptdt(nadv, dtbon, n, pbaru, masse)
Definition: adaptdt.F:6
llmp1
!$Header llmp1
Definition: paramet.h:14
infotrac::ok_iso_verif
logical, save ok_iso_verif
Definition: infotrac.F90:44
llm
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm day day1 day day1 *dt_dice endif!time annee_ref dt_dice swup_dice vg_dice omega_dice tg_prof vg_profd w_profd omega_profd!do llm!print llm l llm
Definition: 1D_interp_cases.h:103
advn
subroutine advn(q, masse, w, pbaru, pbarv, pdt, mode)
Definition: advn.F:5
control_mod::day_step
integer, save day_step
Definition: control_mod.F90:15
infotrac::nqtot
integer, save nqtot
Definition: infotrac.F90:6
scopy
subroutine scopy(n, sx, incx, sy, incy)
Definition: cray.F:9
ip1jm
!$Header llmm1 INTEGER ip1jm
Definition: paramet.h:14
massbar
subroutine massbar(masse, massebx, masseby)
Definition: massbar.F90:2
groupe
subroutine groupe(pext, pbaru, pbarv, pbarum, pbarvm, wm)
Definition: groupe.F:5
jjp1
!$Header jjp1
Definition: paramet.h:14
minmax
subroutine minmax(imax, xi, zmin, zmax)
Definition: minmax.F:5
ppm3d
subroutine ppm3d(IGD, Q, PS1, PS2, U, V, W, NDT, IORD, JORD, KORD, NC, IMR, JNP, j1, NLAY, AP, BP, PT, AE, fill, dum, Umax)
Definition: ppm3d.F:67
prather
subroutine prather(q, w, masse, pbaru, pbarv, nt, dt)
Definition: prather.F:5
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
vlspltqs
subroutine vlspltqs(q, pente_max, masse, w, pbaru, pbarv, pdt, p, pk, teta, iq)
Definition: vlspltqs.F:6
dtvr
!$Id mode_top_bound COMMON comconstr dtvr
Definition: comconst.h:7
check_isotopes_seq
subroutine check_isotopes_seq(q, ip1jmp1, err_msg)
Definition: check_isotopes.F:2
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
infotrac::nqperes
integer, save nqperes
Definition: infotrac.F90:15
vlsplt
subroutine vlsplt(q, pente_max, masse, w, pbaru, pbarv, pdt, iq)
Definition: vlsplt.F:6
advtrac
subroutine advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, flxw, pk)
Definition: advtrac.F90:4
control_mod
Definition: control_mod.F90:5
infotrac::iadv
integer, dimension(:), allocatable, save iadv
Definition: infotrac.F90:22
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7