Context Navigation

advtrac.F90 @ 5160

Last change on this file since 5160 was 1549, checked in by lguez, 13 years ago
Conversion to free source form, no other change
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.1 KB

Line
1	! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z abarral $
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
12	USE control_mod
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	! Ajout PPM
33	!--------------------------------------------------------
34	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
35	!--------------------------------------------------------
36	INTEGER iapptrac
37	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
38	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
39	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
40	REAL pk(ip1jmp1,llm)
41	REAL flxw(ip1jmp1,llm)
42
43	!-------------------------------------------------------------
44	! Variables locales
45	!-------------------------------------------------------------
46
47	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
48	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
49	REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
50	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
51	INTEGER iadvtr
52	INTEGER ij,l,iq,iiq
53	REAL zdpmin, zdpmax
54	EXTERNAL minmax
55	SAVE iadvtr, massem, pbaruc, pbarvc
56	DATA iadvtr/0/
57	!----------------------------------------------------------
58	! Rajouts pour PPM
59	!----------------------------------------------------------
60	INTEGER indice,n
61	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
62	REAL CFLmaxz,aaa,bbb ! CFL maximum
63	REAL psppm(iim,jjp1) ! pression au sol
64	REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
65	REAL qppm(iim*jjp1,llm,nqtot)
66	REAL fluxwppm(iim,jjp1,llm)
67	REAL apppm(llmp1), bpppm(llmp1)
68	LOGICAL dum,fill
69	DATA fill/.true./
70	DATA dum/.true./
71
72	integer,save :: countcfl=0
73	real cflx(ip1jmp1,llm)
74	real cfly(ip1jm,llm)
75	real cflz(ip1jmp1,llm)
76	real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
77
78	IF(iadvtr.EQ.0) THEN
79	CALL initial0(ijp1llm,pbaruc)
80	CALL initial0(ijmllm,pbarvc)
81	ENDIF
82
83	! accumulation des flux de masse horizontaux
84	DO l=1,llm
85	DO ij = 1,ip1jmp1
86	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
87	ENDDO
88	DO ij = 1,ip1jm
89	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
90	ENDDO
91	ENDDO
92
93	! selection de la masse instantannee des mailles avant le transport.
94	IF(iadvtr.EQ.0) THEN
95
96	CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
97	!cc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
98	!
99	ENDIF
100
101	iadvtr = iadvtr+1
102	iapptrac = iadvtr
103
104
105	! Test pour savoir si on advecte a ce pas de temps
106	IF ( iadvtr.EQ.iapp_tracvl ) THEN
107
108	!c .. Modif P.Le Van ( 20/12/97 ) ....
109	!c
110
111	! traitement des flux de masse avant advection.
112	! 1. calcul de w
113	! 2. groupement des mailles pres du pole.
114
115	CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
116
117	! ... Flux de masse diaganostiques traceurs
118	flxw = wg / REAL(iapp_tracvl)
119
120	! test sur l'eventuelle creation de valeurs negatives de la masse
121	DO l=1,llm-1
122	DO ij = iip2+1,ip1jm
123	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l) &
124	- pbarvg(ij-iip1,l) + pbarvg(ij,l) &
125	+ wg(ij,l+1) - wg(ij,l)
126	ENDDO
127	CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
128	DO ij = iip2,ip1jm
129	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
130	ENDDO
131
132
133	CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
134
135	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
136	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin, &
137	' MAX:', zdpmax
138	ENDIF
139
140	ENDDO
141
142
143	!-------------------------------------------------------------------
144	! Calcul des criteres CFL en X, Y et Z
145	!-------------------------------------------------------------------
146
147	if (countcfl == 0. ) then
148	cflxmax(:)=0.
149	cflymax(:)=0.
150	cflzmax(:)=0.
151	endif
152
153	countcfl=countcfl+iapp_tracvl
154	cflx(:,:)=0.
155	cfly(:,:)=0.
156	cflz(:,:)=0.
157	do l=1,llm
158	do ij=iip2,ip1jm-1
159	if (pbarug(ij,l)>=0.) then
160	cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
161	else
162	cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
163	endif
164	enddo
165	enddo
166	do l=1,llm
167	do ij=iip2,ip1jm-1,iip1
168	cflx(ij+iip1,l)=cflx(ij,l)
169	enddo
170	enddo
171
172	do l=1,llm
173	do ij=1,ip1jm
174	if (pbarvg(ij,l)>=0.) then
175	cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
176	else
177	cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
178	endif
179	enddo
180	enddo
181
182	do l=2,llm
183	do ij=1,ip1jm
184	if (wg(ij,l)>=0.) then
185	cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
186	else
187	cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
188	endif
189	enddo
190	enddo
191
192	do l=1,llm
193	cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
194	cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
195	cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
196	enddo
197
198	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
199	! Par defaut, on sort le diagnostic des CFL tous les jours.
200	! Si on veut le sortir a chaque pas d'advection en cas de plantage
201	! if (countcfl==iapp_tracvl) then
202	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203	if (countcfl==day_step) then
204	do l=1,llm
205	write(lunout,*) 'L, CFLmax ' &
206	,l,maxval(cflx(:,l)),maxval(cfly(:,l)),maxval(cflz(:,l))
207	enddo
208	countcfl=0
209	endif
210
211	!-------------------------------------------------------------------
212	! Advection proprement dite (Modification Le Croller (07/2001)
213	!-------------------------------------------------------------------
214
215	!----------------------------------------------------
216	! Calcul des moyennes basées sur la masse
217	!----------------------------------------------------
218	call massbar(massem,massebx,masseby)
219
220	!-----------------------------------------------------------
221	! Appel des sous programmes d'advection
222	!-----------------------------------------------------------
223	do iq=1,nqtot
224	! call clock(t_initial)
225	if(iadv(iq) == 0) cycle
226	! ----------------------------------------------------------------
227	! Schema de Van Leer I MUSCL
228	! ----------------------------------------------------------------
229	if(iadv(iq).eq.10) THEN
230	call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
231
232	! ----------------------------------------------------------------
233	! Schema "pseudo amont" + test sur humidite specifique
234	! pour la vapeur d'eau. F. Codron
235	! ----------------------------------------------------------------
236	else if(iadv(iq).eq.14) then
237	!
238	CALL vlspltqs( q(1,1,1), 2., massem, wg , &
239	pbarug,pbarvg,dtvr,p,pk,teta )
240	! ----------------------------------------------------------------
241	! Schema de Frederic Hourdin
242	! ----------------------------------------------------------------
243	else if(iadv(iq).eq.12) then
244	! Pas de temps adaptatif
245	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
246	if (n.GT.1) then
247	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
248	dtvr,'n=',n
249	endif
250	do indice=1,n
251	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
252	end do
253	else if(iadv(iq).eq.13) then
254	! Pas de temps adaptatif
255	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
256	if (n.GT.1) then
257	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
258	dtvr,'n=',n
259	endif
260	do indice=1,n
261	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
262	end do
263	! ----------------------------------------------------------------
264	! Schema de pente SLOPES
265	! ----------------------------------------------------------------
266	else if (iadv(iq).eq.20) then
267	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
268
269	! ----------------------------------------------------------------
270	! Schema de Prather
271	! ----------------------------------------------------------------
272	else if (iadv(iq).eq.30) then
273	! Pas de temps adaptatif
274	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
275	if (n.GT.1) then
276	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
277	dtvr,'n=',n
278	endif
279	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
280	n,dtbon)
281
282	! ----------------------------------------------------------------
283	! Schemas PPM Lin et Rood
284	! ----------------------------------------------------------------
285	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
286	iadv(iq).LE.18)) then
287
288	! Test sur le flux horizontal
289	! Pas de temps adaptatif
290	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
291	if (n.GT.1) then
292	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=', &
293	dtvr,'n=',n
294	endif
295	! Test sur le flux vertical
296	CFLmaxz=0.
297	do l=2,llm
298	do ij=iip2,ip1jm
299	aaa=wg(ij,l)*dtvr/massem(ij,l)
300	CFLmaxz=max(CFLmaxz,aaa)
301	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
302	CFLmaxz=max(CFLmaxz,bbb)
303	enddo
304	enddo
305	if (CFLmaxz.GE.1) then
306	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
307	endif
308
309	!-----------------------------------------------------------
310	! Ss-prg interface LMDZ.4->PPM3d
311	!-----------------------------------------------------------
312
313	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
314	apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
315	unatppm,vnatppm,psppm)
316
317	do indice=1,n
318	!----------------------------------------------------------------
319	! VL (version PPM) horiz. et PPM vert.
320	!----------------------------------------------------------------
321	if (iadv(iq).eq.11) then
322	! Ss-prg PPM3d de Lin
323	call ppm3d(1,qppm(1,1,iq), &
324	psppm,psppm, &
325	unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
326	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
327	fill,dum,220.)
328
329	!-------------------------------------------------------------
330	! Monotonic PPM
331	!-------------------------------------------------------------
332	else if (iadv(iq).eq.16) then
333	! Ss-prg PPM3d de Lin
334	call ppm3d(1,qppm(1,1,iq), &
335	psppm,psppm, &
336	unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
337	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
338	fill,dum,220.)
339	!-------------------------------------------------------------
340
341	!-------------------------------------------------------------
342	! Semi Monotonic PPM
343	!-------------------------------------------------------------
344	else if (iadv(iq).eq.17) then
345	! Ss-prg PPM3d de Lin
346	call ppm3d(1,qppm(1,1,iq), &
347	psppm,psppm, &
348	unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
349	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
350	fill,dum,220.)
351	!-------------------------------------------------------------
352
353	!-------------------------------------------------------------
354	! Positive Definite PPM
355	!-------------------------------------------------------------
356	else if (iadv(iq).eq.18) then
357	! Ss-prg PPM3d de Lin
358	call ppm3d(1,qppm(1,1,iq), &
359	psppm,psppm, &
360	unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
361	iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
362	fill,dum,220.)
363	!-------------------------------------------------------------
364	endif
365	enddo
366	!-----------------------------------------------------------------
367	! Ss-prg interface PPM3d-LMDZ.4
368	!-----------------------------------------------------------------
369	call interpost(q(1,1,iq),qppm(1,1,iq))
370	endif
371	!----------------------------------------------------------------------
372
373	!-----------------------------------------------------------------
374	! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
375	! et Nord j=1
376	!-----------------------------------------------------------------
377
378	! call traceurpole(q(1,1,iq),massem)
379
380	! calcul du temps cpu pour un schema donne
381
382	! call clock(t_final)
383	!ym tps_cpu=t_final-t_initial
384	!ym cpuadv(iq)=cpuadv(iq)+tps_cpu
385
386	end DO
387
388
389	!------------------------------------------------------------------
390	! on reinitialise a zero les flux de masse cumules
391	!---------------------------------------------------
392	iadvtr=0
393
394	ENDIF ! if iadvtr.EQ.iapp_tracvl
395
396	END SUBROUTINE advtrac

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

Download in other formats:

Original Format