Context Navigation

advtrac_p.F @ 2678

Last change on this file since 2678 was 1114, checked in by jghattas, 16 years ago

Creation du module infotrac:

contient les variables de advtrac.h
contient la subroutine iniadvtrac renommer en infotrac_init
le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
ce module est utilise dans la dynamique et la physique
contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 16.5 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c
6	SUBROUTINE advtrac_p(pbaru,pbarv ,
7	* p, masse,q,iapptrac,teta,
8	* flxw,
9	* pk )
10
11	c Auteur : F. Hourdin
12	c
13	c Modif. P. Le Van (20/12/97)
14	c F. Codron (10/99)
15	c D. Le Croller (07/2001)
16	c M.A Filiberti (04/2002)
17	c
18	USE parallel
19	USE Write_Field_p
20	USE Bands
21	USE mod_hallo
22	USE Vampir
23	USE times
24	USE infotrac
25	IMPLICIT NONE
26	c
27	#include "dimensions.h"
28	#include "paramet.h"
29	#include "comconst.h"
30	#include "comvert.h"
31	#include "comdissip.h"
32	#include "comgeom2.h"
33	#include "logic.h"
34	#include "temps.h"
35	#include "control.h"
36	#include "ener.h"
37	#include "description.h"
38
39	c-------------------------------------------------------------------
40	c Arguments
41	c-------------------------------------------------------------------
42	c Ajout PPM
43	c--------------------------------------------------------
44	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
45	c--------------------------------------------------------
46	INTEGER iapptrac
47	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
48	REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
49	REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
50	REAL pk(ip1jmp1,llm)
51	REAL :: flxw(ip1jmp1,llm)
52
53	c-------------------------------------------------------------
54	c Variables locales
55	c-------------------------------------------------------------
56
57	REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
58	REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
59	REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
60	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
61	INTEGER iadvtr
62	INTEGER ij,l,iq,iiq
63	REAL zdpmin, zdpmax
64	SAVE iadvtr, massem, pbaruc, pbarvc
65	DATA iadvtr/0/
66	c$OMP THREADPRIVATE(iadvtr)
67	c----------------------------------------------------------
68	c Rajouts pour PPM
69	c----------------------------------------------------------
70	INTEGER indice,n
71	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
72	REAL CFLmaxz,aaa,bbb ! CFL maximum
73	REAL psppm(iim,jjp1) ! pression au sol
74	REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
75	REAL qppm(iim*jjp1,llm,nqtot)
76	REAL fluxwppm(iim,jjp1,llm)
77	REAL apppm(llmp1), bpppm(llmp1)
78	LOGICAL dum,fill
79	DATA fill/.true./
80	DATA dum/.true./
81	REAL,SAVE :: finmasse(ip1jmp1,llm)
82	integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
83	type(Request) :: Request_vanleer
84	REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
85	REAL,SAVE :: teta_tmp(ip1jmp1,llm)
86	REAL,SAVE :: pk_tmp(ip1jmp1,llm)
87
88	ijb_u=ij_begin
89	ije_u=ij_end
90
91	ijb_v=ij_begin-iip1
92	ije_v=ij_end
93	if (pole_nord) ijb_v=ij_begin
94	if (pole_sud) ije_v=ij_end-iip1
95
96	IF(iadvtr.EQ.0) THEN
97	c CALL initial0(ijp1llm,pbaruc)
98	c CALL initial0(ijmllm,pbarvc)
99	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
100	DO l=1,llm
101	pbaruc(ijb_u:ije_u,l)=0.
102	pbarvc(ijb_v:ije_v,l)=0.
103	ENDDO
104	c$OMP END DO NOWAIT
105	ENDIF
106
107	c accumulation des flux de masse horizontaux
108	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
109	DO l=1,llm
110	DO ij = ijb_u,ije_u
111	pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
112	ENDDO
113	DO ij = ijb_v,ije_v
114	pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
115	ENDDO
116	ENDDO
117	c$OMP END DO NOWAIT
118
119	c selection de la masse instantannee des mailles avant le transport.
120	IF(iadvtr.EQ.0) THEN
121
122	c CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
123	ijb=ij_begin
124	ije=ij_end
125
126	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
127	DO l=1,llm
128	massem(ijb:ije,l)=masse(ijb:ije,l)
129	ENDDO
130	c$OMP END DO NOWAIT
131
132	ccc CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
133	c
134	ENDIF
135
136	iadvtr = iadvtr+1
137
138	c$OMP MASTER
139	iapptrac = iadvtr
140	c$OMP END MASTER
141
142	c Test pour savoir si on advecte a ce pas de temps
143
144	IF ( iadvtr.EQ.iapp_tracvl ) THEN
145	c$OMP MASTER
146	call suspend_timer(timer_caldyn)
147	c$OMP END MASTER
148
149	ijb=ij_begin
150	ije=ij_end
151
152
153	cc .. Modif P.Le Van ( 20/12/97 ) ....
154	cc
155
156	c traitement des flux de masse avant advection.
157	c 1. calcul de w
158	c 2. groupement des mailles pres du pole.
159
160	CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
161
162	c$OMP BARRIER
163
164	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
165	DO l=1,llmp1
166	p_tmp(ijb:ije,l)=p(ijb:ije,l)
167	ENDDO
168	c$OMP END DO NOWAIT
169
170	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
171	DO l=1,llm
172	pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
173	teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
174	ENDDO
175	c$OMP END DO NOWAIT
176
177	c$OMP MASTER
178	call VTb(VTHallo)
179	c$OMP END MASTER
180
181	call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
182	* jj_Nb_vanleer,0,0,Request_vanleer)
183	call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
184	* jj_Nb_vanleer,1,0,Request_vanleer)
185	call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
186	* jj_Nb_vanleer,0,0,Request_vanleer)
187	call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
188	* jj_Nb_vanleer,0,0,Request_vanleer)
189	call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
190	* jj_Nb_vanleer,1,1,Request_vanleer)
191	call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
192	* jj_Nb_vanleer,1,1,Request_vanleer)
193	call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
194	* jj_Nb_vanleer,1,1,Request_vanleer)
195	do j=1,nqtot
196	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
197	* jj_nb_vanleer,0,0,Request_vanleer)
198	enddo
199
200	call SendRequest(Request_vanleer)
201	c$OMP BARRIER
202	call WaitRequest(Request_vanleer)
203
204
205	c$OMP BARRIER
206	c$OMP MASTER
207	call SetDistrib(jj_nb_vanleer)
208	call VTe(VTHallo)
209	call VTb(VTadvection)
210	call start_timer(timer_vanleer)
211	c$OMP END MASTER
212	c$OMP BARRIER
213
214	! ... Flux de masse diaganostiques traceurs
215	ijb=ij_begin
216	ije=ij_end
217	flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/FLOAT(iapp_tracvl)
218
219	c test sur l'eventuelle creation de valeurs negatives de la masse
220	ijb=ij_begin
221	ije=ij_end
222	if (pole_nord) ijb=ij_begin+iip1
223	if (pole_sud) ije=ij_end-iip1
224
225	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
226	DO l=1,llm-1
227	DO ij = ijb+1,ije
228	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
229	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
230	s + wg(ij,l+1) - wg(ij,l)
231	ENDDO
232
233	c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
234	c ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
235
236	do ij=ijb,ije-iip1+1,iip1
237	zdp(ij)=zdp(ij+iip1-1)
238	enddo
239
240	DO ij = ijb,ije
241	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
242	ENDDO
243
244
245	c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
246	c ym ---> eventuellement a revoir
247	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
248
249	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
250	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
251	s ' MAX:', zdpmax
252	ENDIF
253
254	ENDDO
255	c$OMP END DO NOWAIT
256
257	c-------------------------------------------------------------------
258	c Advection proprement dite (Modification Le Croller (07/2001)
259	c-------------------------------------------------------------------
260
261	c----------------------------------------------------
262	c Calcul des moyennes basées sur la masse
263	c----------------------------------------------------
264
265	cym ----> Normalement, inutile pour les schémas classiques
266	cym ----> Revérifier lors de la parallélisation des autres schemas
267
268	cym call massbar_p(massem,massebx,masseby)
269
270	call vlspltgen_p( q,iadv, 2., massem, wg ,
271	* pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
272
273
274	GOTO 1234
275	c-----------------------------------------------------------
276	c Appel des sous programmes d'advection
277	c-----------------------------------------------------------
278	do iq=1,nqtot
279	c call clock(t_initial)
280	if(iadv(iq) == 0) cycle
281	c ----------------------------------------------------------------
282	c Schema de Van Leer I MUSCL
283	c ----------------------------------------------------------------
284	if(iadv(iq).eq.10) THEN
285
286	call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
287
288	c ----------------------------------------------------------------
289	c Schema "pseudo amont" + test sur humidite specifique
290	C pour la vapeur d'eau. F. Codron
291	c ----------------------------------------------------------------
292	else if(iadv(iq).eq.14) then
293	c
294	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
295	CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
296	* pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
297	c ----------------------------------------------------------------
298	c Schema de Frederic Hourdin
299	c ----------------------------------------------------------------
300	else if(iadv(iq).eq.12) then
301	stop 'advtrac : schema non parallelise'
302	c 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	s dtvr,'n=',n
307	endif
308	do indice=1,n
309	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
310	end do
311	else if(iadv(iq).eq.13) then
312	stop 'advtrac : schema non parallelise'
313	c Pas de temps adaptatif
314	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
315	if (n.GT.1) then
316	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
317	s dtvr,'n=',n
318	endif
319	do indice=1,n
320	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
321	end do
322	c ----------------------------------------------------------------
323	c Schema de pente SLOPES
324	c ----------------------------------------------------------------
325	else if (iadv(iq).eq.20) then
326	stop 'advtrac : schema non parallelise'
327
328	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
329
330	c ----------------------------------------------------------------
331	c Schema de Prather
332	c ----------------------------------------------------------------
333	else if (iadv(iq).eq.30) then
334	stop 'advtrac : schema non parallelise'
335	c Pas de temps adaptatif
336	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
337	if (n.GT.1) then
338	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
339	s dtvr,'n=',n
340	endif
341	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
342	s n,dtbon)
343	c ----------------------------------------------------------------
344	c Schemas PPM Lin et Rood
345	c ----------------------------------------------------------------
346	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
347	s iadv(iq).LE.18)) then
348
349	stop 'advtrac : schema non parallelise'
350
351	c Test sur le flux horizontal
352	c Pas de temps adaptatif
353	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
354	if (n.GT.1) then
355	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
356	s dtvr,'n=',n
357	endif
358	c Test sur le flux vertical
359	CFLmaxz=0.
360	do l=2,llm
361	do ij=iip2,ip1jm
362	aaa=wg(ij,l)*dtvr/massem(ij,l)
363	CFLmaxz=max(CFLmaxz,aaa)
364	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
365	CFLmaxz=max(CFLmaxz,bbb)
366	enddo
367	enddo
368	if (CFLmaxz.GE.1) then
369	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
370	endif
371
372	c-----------------------------------------------------------
373	c Ss-prg interface LMDZ.4->PPM3d
374	c-----------------------------------------------------------
375
376	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
377	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
378	s unatppm,vnatppm,psppm)
379
380	do indice=1,n
381	c---------------------------------------------------------------------
382	c VL (version PPM) horiz. et PPM vert.
383	c---------------------------------------------------------------------
384	if (iadv(iq).eq.11) then
385	c Ss-prg PPM3d de Lin
386	call ppm3d(1,qppm(1,1,iq),
387	s psppm,psppm,
388	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
389	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
390	s fill,dum,220.)
391
392	c----------------------------------------------------------------------
393	c Monotonic PPM
394	c----------------------------------------------------------------------
395	else if (iadv(iq).eq.16) then
396	c Ss-prg PPM3d de Lin
397	call ppm3d(1,qppm(1,1,iq),
398	s psppm,psppm,
399	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
400	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
401	s fill,dum,220.)
402	c---------------------------------------------------------------------
403
404	c---------------------------------------------------------------------
405	c Semi Monotonic PPM
406	c---------------------------------------------------------------------
407	else if (iadv(iq).eq.17) then
408	c Ss-prg PPM3d de Lin
409	call ppm3d(1,qppm(1,1,iq),
410	s psppm,psppm,
411	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
412	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
413	s fill,dum,220.)
414	c---------------------------------------------------------------------
415
416	c---------------------------------------------------------------------
417	c Positive Definite PPM
418	c---------------------------------------------------------------------
419	else if (iadv(iq).eq.18) then
420	c Ss-prg PPM3d de Lin
421	call ppm3d(1,qppm(1,1,iq),
422	s psppm,psppm,
423	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
424	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
425	s fill,dum,220.)
426	c---------------------------------------------------------------------
427	endif
428	enddo
429	c-----------------------------------------------------------------
430	c Ss-prg interface PPM3d-LMDZ.4
431	c-----------------------------------------------------------------
432	call interpost(q(1,1,iq),qppm(1,1,iq))
433	endif
434	c----------------------------------------------------------------------
435
436	c-----------------------------------------------------------------
437	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
438	c et Nord j=1
439	c-----------------------------------------------------------------
440
441	c call traceurpole(q(1,1,iq),massem)
442
443	c calcul du temps cpu pour un schema donne
444
445	c call clock(t_final)
446	cym tps_cpu=t_final-t_initial
447	cym cpuadv(iq)=cpuadv(iq)+tps_cpu
448
449	end DO
450
451	1234 CONTINUE
452	c$OMP BARRIER
453
454	ijb=ij_begin
455	ije=ij_end
456
457	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
458	DO l = 1, llm
459	DO ij = ijb, ije
460	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
461	ENDDO
462	ENDDO
463	c$OMP END DO
464
465	CALL qminimum_p( q, 2, finmasse )
466
467	c------------------------------------------------------------------
468	c on reinitialise a zero les flux de masse cumules
469	c---------------------------------------------------
470	c iadvtr=0
471
472	c$OMP MASTER
473	call VTe(VTadvection)
474	call stop_timer(timer_vanleer)
475	call VTb(VThallo)
476	c$OMP END MASTER
477
478	do j=1,nqtot
479	call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
480	* jj_nb_caldyn,0,0,Request_vanleer)
481	enddo
482
483	call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
484	* jj_nb_caldyn,0,0,Request_vanleer)
485
486	call SendRequest(Request_vanleer)
487	c$OMP BARRIER
488	call WaitRequest(Request_vanleer)
489
490	c$OMP BARRIER
491	c$OMP MASTER
492	call SetDistrib(jj_nb_caldyn)
493	call VTe(VThallo)
494	call resume_timer(timer_caldyn)
495	c$OMP END MASTER
496	c$OMP BARRIER
497	iadvtr=0
498	ENDIF ! if iadvtr.EQ.iapp_tracvl
499
500	RETURN
501	END
502

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

Download in other formats:

Original Format