Context Navigation

advtrac.F @ 5000

Last change on this file since 5000 was 1176, checked in by Laurent Fairhead, 15 years ago
Modif pour la compilation avec fcm sur Vargas SD Corrections relatives au controle de histins et de l'utilisation des options INST(X) pour l'appel au routines de IOIPSL FH Mise en place de diagnostics sur les critères CFL pour l'advection de traceurs. FH Controle dans les .def de la dépendance verticale de l'efficacite de la diffusion. Actif pour le moment uniquement avec ok_strato=y et llm=39. FH LF
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 1176 2009-06-11 08:54:10Z abarral $
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