Context Navigation

advtrac.F90 @ 4196

Last change on this file since 4196 was 2286, checked in by Ehouarn Millour, 9 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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 2286 2015-05-20 13:27:07Z fairhead $
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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format