Context Navigation

advtrac.F @ 3226

Last change on this file since 3226 was 1454, checked in by Laurent Fairhead, 14 years ago

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434

Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

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 1454 2010-11-18 12:01:24Z fairhead $
3	!
4	c
5	c
6	SUBROUTINE advtrac(pbaru,pbarv ,
7	* p, masse,q,iapptrac,teta,
8	* flxw,
9	* pk)
10	c Auteur : F. Hourdin
11	c
12	c Modif. P. Le Van (20/12/97)
13	c F. Codron (10/99)
14	c D. Le Croller (07/2001)
15	c M.A Filiberti (04/2002)
16	c
17	USE infotrac
18	USE control_mod
19
20
21	IMPLICIT NONE
22	c
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
35	c-------------------------------------------------------------------
36	c Arguments
37	c-------------------------------------------------------------------
38	c Ajout PPM
39	c--------------------------------------------------------
40	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
41	c--------------------------------------------------------
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
49	c-------------------------------------------------------------
50	c Variables locales
51	c-------------------------------------------------------------
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/
63	c----------------------------------------------------------
64	c Rajouts pour PPM
65	c----------------------------------------------------------
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
89	c 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
99	c 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)
103	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
104	c
105	ENDIF
106
107	iadvtr = iadvtr+1
108	iapptrac = iadvtr
109
110
111	c Test pour savoir si on advecte a ce pas de temps
112	IF ( iadvtr.EQ.iapp_tracvl ) THEN
113
114	cc .. Modif P.Le Van ( 20/12/97 ) ....
115	cc
116
117	c traitement des flux de masse avant advection.
118	c 1. calcul de w
119	c 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
126	c 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
149	c-------------------------------------------------------------------
150	! Calcul des criteres CFL en X, Y et Z
151	c-------------------------------------------------------------------
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
217	c-------------------------------------------------------------------
218	c Advection proprement dite (Modification Le Croller (07/2001)
219	c-------------------------------------------------------------------
220
221	c----------------------------------------------------
222	c Calcul des moyennes basées sur la masse
223	c----------------------------------------------------
224	call massbar(massem,massebx,masseby)
225
226	c-----------------------------------------------------------
227	c Appel des sous programmes d'advection
228	c-----------------------------------------------------------
229	do iq=1,nqtot
230	c call clock(t_initial)
231	if(iadv(iq) == 0) cycle
232	c ----------------------------------------------------------------
233	c Schema de Van Leer I MUSCL
234	c ----------------------------------------------------------------
235	if(iadv(iq).eq.10) THEN
236	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
237
238	c ----------------------------------------------------------------
239	c Schema "pseudo amont" + test sur humidite specifique
240	C pour la vapeur d'eau. F. Codron
241	c ----------------------------------------------------------------
242	else if(iadv(iq).eq.14) then
243	c
244	CALL vlspltqs( q(1,1,1), 2., massem, wg ,
245	* pbarug,pbarvg,dtvr,p,pk,teta )
246	c ----------------------------------------------------------------
247	c Schema de Frederic Hourdin
248	c ----------------------------------------------------------------
249	else if(iadv(iq).eq.12) then
250	c Pas de temps adaptatif
251	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
252	if (n.GT.1) then
253	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
254	s dtvr,'n=',n
255	endif
256	do indice=1,n
257	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
258	end do
259	else if(iadv(iq).eq.13) then
260	c Pas de temps adaptatif
261	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
262	if (n.GT.1) then
263	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
264	s dtvr,'n=',n
265	endif
266	do indice=1,n
267	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
268	end do
269	c ----------------------------------------------------------------
270	c Schema de pente SLOPES
271	c ----------------------------------------------------------------
272	else if (iadv(iq).eq.20) then
273	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
274
275	c ----------------------------------------------------------------
276	c Schema de Prather
277	c ----------------------------------------------------------------
278	else if (iadv(iq).eq.30) then
279	c Pas de temps adaptatif
280	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
281	if (n.GT.1) then
282	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
283	s dtvr,'n=',n
284	endif
285	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
286	s n,dtbon)
287
288	c ----------------------------------------------------------------
289	c Schemas PPM Lin et Rood
290	c ----------------------------------------------------------------
291	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
292	s iadv(iq).LE.18)) then
293
294	c Test sur le flux horizontal
295	c Pas de temps adaptatif
296	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
297	if (n.GT.1) then
298	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
299	s dtvr,'n=',n
300	endif
301	c Test sur le flux vertical
302	CFLmaxz=0.
303	do l=2,llm
304	do ij=iip2,ip1jm
305	aaa=wg(ij,l)*dtvr/massem(ij,l)
306	CFLmaxz=max(CFLmaxz,aaa)
307	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
308	CFLmaxz=max(CFLmaxz,bbb)
309	enddo
310	enddo
311	if (CFLmaxz.GE.1) then
312	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
313	endif
314
315	c-----------------------------------------------------------
316	c Ss-prg interface LMDZ.4->PPM3d
317	c-----------------------------------------------------------
318
319	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
320	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
321	s unatppm,vnatppm,psppm)
322
323	do indice=1,n
324	c---------------------------------------------------------------------
325	c VL (version PPM) horiz. et PPM vert.
326	c---------------------------------------------------------------------
327	if (iadv(iq).eq.11) then
328	c Ss-prg PPM3d de Lin
329	call ppm3d(1,qppm(1,1,iq),
330	s psppm,psppm,
331	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
332	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
333	s fill,dum,220.)
334
335	c----------------------------------------------------------------------
336	c Monotonic PPM
337	c----------------------------------------------------------------------
338	else if (iadv(iq).eq.16) then
339	c Ss-prg PPM3d de Lin
340	call ppm3d(1,qppm(1,1,iq),
341	s psppm,psppm,
342	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
343	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
344	s fill,dum,220.)
345	c---------------------------------------------------------------------
346
347	c---------------------------------------------------------------------
348	c Semi Monotonic PPM
349	c---------------------------------------------------------------------
350	else if (iadv(iq).eq.17) then
351	c Ss-prg PPM3d de Lin
352	call ppm3d(1,qppm(1,1,iq),
353	s psppm,psppm,
354	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
355	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
356	s fill,dum,220.)
357	c---------------------------------------------------------------------
358
359	c---------------------------------------------------------------------
360	c Positive Definite PPM
361	c---------------------------------------------------------------------
362	else if (iadv(iq).eq.18) then
363	c Ss-prg PPM3d de Lin
364	call ppm3d(1,qppm(1,1,iq),
365	s psppm,psppm,
366	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
367	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
368	s fill,dum,220.)
369	c---------------------------------------------------------------------
370	endif
371	enddo
372	c-----------------------------------------------------------------
373	c Ss-prg interface PPM3d-LMDZ.4
374	c-----------------------------------------------------------------
375	call interpost(q(1,1,iq),qppm(1,1,iq))
376	endif
377	c----------------------------------------------------------------------
378
379	c-----------------------------------------------------------------
380	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
381	c et Nord j=1
382	c-----------------------------------------------------------------
383
384	c call traceurpole(q(1,1,iq),massem)
385
386	c calcul du temps cpu pour un schema donne
387
388	c call clock(t_final)
389	cym tps_cpu=t_final-t_initial
390	cym cpuadv(iq)=cpuadv(iq)+tps_cpu
391
392	end DO
393
394
395	c------------------------------------------------------------------
396	c on reinitialise a zero les flux de masse cumules
397	c---------------------------------------------------
398	iadvtr=0
399
400	ENDIF ! if iadvtr.EQ.iapp_tracvl
401
402	RETURN
403	END
404

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

Download in other formats:

Original Format