Context Navigation

advtrac.F

Last change on this file was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.7 KB

Line
1	!
2	! $Id: advtrac.F 1279 2009-12-10 09:02:56Z aclsce $
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
19	IMPLICIT NONE
20	c
21	#include "dimensions.h"
22	#include "paramet.h"
23	#include "comconst.h"
24	#include "comvert.h"
25	#include "comdissip.h"
26	#include "comgeom2.h"
27	#include "logic.h"
28	#include "temps.h"
29	#include "control.h"
30	#include "ener.h"
31	#include "description.h"
32	#include "iniprint.h"
33
34	c-------------------------------------------------------------------
35	c Arguments
36	c-------------------------------------------------------------------
37	c Ajout PPM
38	c--------------------------------------------------------
39	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
40	c--------------------------------------------------------
41	INTEGER iapptrac
42	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
43	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
44	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
45	REAL pk(ip1jmp1,llm)
46	REAL flxw(ip1jmp1,llm)
47
48	c-------------------------------------------------------------
49	c Variables locales
50	c-------------------------------------------------------------
51
52	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
53	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
54	REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
55	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
56	INTEGER iadvtr
57	INTEGER ij,l,iq,iiq
58	REAL zdpmin, zdpmax
59	EXTERNAL minmax
60	SAVE iadvtr, massem, pbaruc, pbarvc
61	DATA iadvtr/0/
62	c----------------------------------------------------------
63	c Rajouts pour PPM
64	c----------------------------------------------------------
65	INTEGER indice,n
66	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
67	REAL CFLmaxz,aaa,bbb ! CFL maximum
68	REAL psppm(iim,jjp1) ! pression au sol
69	REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
70	REAL qppm(iim*jjp1,llm,nqtot)
71	REAL fluxwppm(iim,jjp1,llm)
72	REAL apppm(llmp1), bpppm(llmp1)
73	LOGICAL dum,fill
74	DATA fill/.true./
75	DATA dum/.true./
76
77	integer,save :: countcfl=0
78	real cflx(ip1jmp1,llm)
79	real cfly(ip1jm,llm)
80	real cflz(ip1jmp1,llm)
81	real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
82
83	IF(iadvtr.EQ.0) THEN
84	CALL initial0(ijp1llm,pbaruc)
85	CALL initial0(ijmllm,pbarvc)
86	ENDIF
87
88	c accumulation des flux de masse horizontaux
89	DO l=1,llm
90	DO ij = 1,ip1jmp1
91	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
92	ENDDO
93	DO ij = 1,ip1jm
94	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
95	ENDDO
96	ENDDO
97
98	c selection de la masse instantannee des mailles avant le transport.
99	IF(iadvtr.EQ.0) THEN
100
101	CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
102	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
103	c
104	ENDIF
105
106	iadvtr = iadvtr+1
107	iapptrac = iadvtr
108
109
110	c Test pour savoir si on advecte a ce pas de temps
111	IF ( iadvtr.EQ.iapp_tracvl ) THEN
112
113	cc .. Modif P.Le Van ( 20/12/97 ) ....
114	cc
115
116	c traitement des flux de masse avant advection.
117	c 1. calcul de w
118	c 2. groupement des mailles pres du pole.
119
120	CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
121
122	! ... Flux de masse diaganostiques traceurs
123	flxw = wg / FLOAT(iapp_tracvl)
124
125	c test sur l'eventuelle creation de valeurs negatives de la masse
126	DO l=1,llm-1
127	DO ij = iip2+1,ip1jm
128	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
129	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
130	s + wg(ij,l+1) - wg(ij,l)
131	ENDDO
132	CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
133	DO ij = iip2,ip1jm
134	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
135	ENDDO
136
137
138	CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
139
140	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
141	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
142	s ' MAX:', zdpmax
143	ENDIF
144
145	ENDDO
146
147
148	c-------------------------------------------------------------------
149	! Calcul des criteres CFL en X, Y et Z
150	c-------------------------------------------------------------------
151
152	if (countcfl == 0. ) then
153	cflxmax(:)=0.
154	cflymax(:)=0.
155	cflzmax(:)=0.
156	endif
157
158	countcfl=countcfl+iapp_tracvl
159	cflx(:,:)=0.
160	cfly(:,:)=0.
161	cflz(:,:)=0.
162	do l=1,llm
163	do ij=iip2,ip1jm-1
164	if (pbarug(ij,l)>=0.) then
165	cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
166	else
167	cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
168	endif
169	enddo
170	enddo
171	do l=1,llm
172	do ij=iip2,ip1jm-1,iip1
173	cflx(ij+iip1,l)=cflx(ij,l)
174	enddo
175	enddo
176
177	do l=1,llm
178	do ij=1,ip1jm
179	if (pbarvg(ij,l)>=0.) then
180	cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
181	else
182	cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
183	endif
184	enddo
185	enddo
186
187	do l=2,llm
188	do ij=1,ip1jm
189	if (wg(ij,l)>=0.) then
190	cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
191	else
192	cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
193	endif
194	enddo
195	enddo
196
197	do l=1,llm
198	cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
199	cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
200	cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
201	enddo
202
203	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
204	! Par defaut, on sort le diagnostic des CFL tous les jours.
205	! Si on veut le sortir a chaque pas d'advection en cas de plantage
206	! if (countcfl==iapp_tracvl) then
207	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208	if (countcfl==day_step) then
209	do l=1,llm
210	write(lunout,*) 'L, CFLmax '
211	s ,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
212	enddo
213	countcfl=0
214	endif
215
216	c-------------------------------------------------------------------
217	c Advection proprement dite (Modification Le Croller (07/2001)
218	c-------------------------------------------------------------------
219
220	c----------------------------------------------------
221	c Calcul des moyennes basées sur la masse
222	c----------------------------------------------------
223	call massbar(massem,massebx,masseby)
224
225	c-----------------------------------------------------------
226	c Appel des sous programmes d'advection
227	c-----------------------------------------------------------
228	do iq=1,nqtot
229	c call clock(t_initial)
230	if(iadv(iq) == 0) cycle
231	c ----------------------------------------------------------------
232	c Schema de Van Leer I MUSCL
233	c ----------------------------------------------------------------
234	if(iadv(iq).eq.10) THEN
235	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
236
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