Context Navigation

advtrac_loc.F @ 4400

Last change on this file since 4400 was 1673, checked in by Laurent Fairhead, 12 years ago

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer

End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

File size: 12.6 KB

Line
1	!
2	! $Id$
3	!
4	c
5	c
6	#define DEBUG_IO
7	#undef DEBUG_IO
8	SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
9	* p, massem,q,teta,
10	* pk )
11
12	c Auteur : F. Hourdin
13	c
14	c Modif. P. Le Van (20/12/97)
15	c F. Codron (10/99)
16	c D. Le Croller (07/2001)
17	c M.A Filiberti (04/2002)
18	c
19	USE parallel
20	USE Write_Field_loc
21	USE Write_Field
22	USE Bands
23	USE mod_hallo
24	USE Vampir
25	USE times
26	USE infotrac
27	USE control_mod
28	USE advtrac_mod
29	IMPLICIT NONE
30	c
31	#include "dimensions.h"
32	#include "paramet.h"
33	#include "comconst.h"
34	#include "comvert.h"
35	#include "comdissip.h"
36	#include "comgeom2.h"
37	#include "logic.h"
38	#include "temps.h"
39	#include "ener.h"
40	#include "description.h"
41
42	c-------------------------------------------------------------------
43	c Arguments
44	c-------------------------------------------------------------------
45	c Ajout PPM
46	c--------------------------------------------------------
47	REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
48	c--------------------------------------------------------
49	INTEGER iapptrac
50	REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
51	REAL wg(ijb_u:ije_u,llm)
52	REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
53	REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
54	REAL pk(ijb_u:ije_u,llm)
55
56	c-------------------------------------------------------------
57	c Variables locales
58	c-------------------------------------------------------------
59
60	REAL zdp(ijb_u:ije_u)
61	REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
62	INTEGER,SAVE :: iadvtr=0
63	c$OMP THREADPRIVATE(iadvtr)
64	INTEGER ij,l,iq,iiq
65	REAL zdpmin, zdpmax
66	c----------------------------------------------------------
67	c Rajouts pour PPM
68	c----------------------------------------------------------
69	INTEGER indice,n
70	REAL dtbon ! Pas de temps adaptatif pour que CFL<1
71	REAL CFLmaxz,aaa,bbb ! CFL maximum
72	REAL psppm(iim,jjb_u:jje_u) ! pression au sol
73	REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
74	REAL qppm(iim*jjnb_u,llm,nqtot)
75	REAL fluxwppm(iim,jjb_u:jje_u,llm)
76	REAL apppm(llmp1), bpppm(llmp1)
77	LOGICAL dum,fill
78	DATA fill/.true./
79	DATA dum/.true./
80	integer ijb,ije,ijbu,ijbv,ijeu,ijev,j
81	type(Request) :: testRequest
82
83	c test sur l'eventuelle creation de valeurs negatives de la masse
84	ijb=ij_begin
85	ije=ij_end
86	if (pole_nord) ijb=ij_begin+iip1
87	if (pole_sud) ije=ij_end-iip1
88
89	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
90	DO l=1,llm-1
91	DO ij = ijb+1,ije
92	zdp(ij) = pbarug(ij-1,l) - pbarug(ij,l)
93	s - pbarvg(ij-iip1,l) + pbarvg(ij,l)
94	s + wg(ij,l+1) - wg(ij,l)
95	ENDDO
96
97	c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
98	c ym ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
99
100	do ij=ijb,ije-iip1+1,iip1
101	zdp(ij)=zdp(ij+iip1-1)
102	enddo
103
104	DO ij = ijb,ije
105	zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
106	ENDDO
107
108
109	c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
110	c ym ---> eventuellement a revoir
111	CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
112
113	IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
114	PRINT*,'WARNING DP/P l=',l,' MIN:',zdpmin,
115	s ' MAX:', zdpmax
116	ENDIF
117
118	ENDDO
119	c$OMP END DO NOWAIT
120
121	c-------------------------------------------------------------------
122	c Advection proprement dite (Modification Le Croller (07/2001)
123	c-------------------------------------------------------------------
124
125	c----------------------------------------------------
126	c Calcul des moyennes basées sur la masse
127	c----------------------------------------------------
128
129	cym ----> Normalement, inutile pour les schémas classiques
130	cym ----> Revérifier lors de la parallélisation des autres schemas
131
132	cym call massbar_p(massem,massebx,masseby)
133
134	#ifdef DEBUG_IO
135	CALL WriteField_u('massem',massem)
136	CALL WriteField_u('wg',wg)
137	CALL WriteField_u('pbarug',pbarug)
138	CALL WriteField_v('pbarvg',pbarvg)
139	CALL WriteField_u('p_tmp',p)
140	CALL WriteField_u('pk_tmp',pk)
141	CALL WriteField_u('teta_tmp',teta)
142	do j=1,nqtot
143	call WriteField_u('q_adv'//trim(int2str(j)),
144	. q(:,:,j))
145	enddo
146	#endif
147
148	!
149	! call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
150	!
151	! call SendRequest(TestRequest)
152	!c$OMP BARRIER
153	! call WaitRequest(TestRequest)
154	c$OMP BARRIER
155
156	call vlspltgen_loc( q,iadv, 2., massem, wg ,
157	* pbarug,pbarvg,dtvr,p,
158	* pk,teta )
159
160	#ifdef DEBUG_IO
161	do j=1,nqtot
162	call WriteField_u('q_adv'//trim(int2str(j)),
163	. q(:,:,j))
164	enddo
165	#endif
166
167	GOTO 1234
168	c-----------------------------------------------------------
169	c Appel des sous programmes d'advection
170	c-----------------------------------------------------------
171	do iq=1,nqtot
172	c call clock(t_initial)
173	if(iadv(iq) == 0) cycle
174	c ----------------------------------------------------------------
175	c Schema de Van Leer I MUSCL
176	c ----------------------------------------------------------------
177	if(iadv(iq).eq.10) THEN
178
179	!LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
180
181	c ----------------------------------------------------------------
182	c Schema "pseudo amont" + test sur humidite specifique
183	C pour la vapeur d'eau. F. Codron
184	c ----------------------------------------------------------------
185	else if(iadv(iq).eq.14) then
186	c
187	cym stop 'advtrac : appel à vlspltqs :schema non parallelise'
188	!LF CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
189	!LF * pbarug,pbarvg,dtvr,p,pk,teta )
190	c ----------------------------------------------------------------
191	c Schema de Frederic Hourdin
192	c ----------------------------------------------------------------
193	else if(iadv(iq).eq.12) then
194	stop 'advtrac : schema non parallelise'
195	c Pas de temps adaptatif
196	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
197	if (n.GT.1) then
198	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
199	s dtvr,'n=',n
200	endif
201	do indice=1,n
202	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
203	end do
204	else if(iadv(iq).eq.13) then
205	stop 'advtrac : schema non parallelise'
206	c Pas de temps adaptatif
207	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
208	if (n.GT.1) then
209	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
210	s dtvr,'n=',n
211	endif
212	do indice=1,n
213	call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
214	end do
215	c ----------------------------------------------------------------
216	c Schema de pente SLOPES
217	c ----------------------------------------------------------------
218	else if (iadv(iq).eq.20) then
219	stop 'advtrac : schema non parallelise'
220
221	call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
222
223	c ----------------------------------------------------------------
224	c Schema de Prather
225	c ----------------------------------------------------------------
226	else if (iadv(iq).eq.30) then
227	stop 'advtrac : schema non parallelise'
228	c Pas de temps adaptatif
229	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
230	if (n.GT.1) then
231	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
232	s dtvr,'n=',n
233	endif
234	call prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
235	s n,dtbon)
236	c ----------------------------------------------------------------
237	c Schemas PPM Lin et Rood
238	c ----------------------------------------------------------------
239	else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
240	s iadv(iq).LE.18)) then
241
242	stop 'advtrac : schema non parallelise'
243
244	c Test sur le flux horizontal
245	c Pas de temps adaptatif
246	call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
247	if (n.GT.1) then
248	write(,) 'WARNING horizontal dt=',dtbon,'dtvr=',
249	s dtvr,'n=',n
250	endif
251	c Test sur le flux vertical
252	CFLmaxz=0.
253	do l=2,llm
254	do ij=iip2,ip1jm
255	aaa=wg(ij,l)*dtvr/massem(ij,l)
256	CFLmaxz=max(CFLmaxz,aaa)
257	bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
258	CFLmaxz=max(CFLmaxz,bbb)
259	enddo
260	enddo
261	if (CFLmaxz.GE.1) then
262	write(,) 'WARNING vertical','CFLmaxz=', CFLmaxz
263	endif
264
265	c-----------------------------------------------------------
266	c Ss-prg interface LMDZ.4->PPM3d
267	c-----------------------------------------------------------
268
269	call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
270	s apppm,bpppm,massebx,masseby,pbarug,pbarvg,
271	s unatppm,vnatppm,psppm)
272
273	do indice=1,n
274	c---------------------------------------------------------------------
275	c VL (version PPM) horiz. et PPM vert.
276	c---------------------------------------------------------------------
277	if (iadv(iq).eq.11) then
278	c Ss-prg PPM3d de Lin
279	call ppm3d(1,qppm(1,1,iq),
280	s psppm,psppm,
281	s unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
282	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
283	s fill,dum,220.)
284
285	c----------------------------------------------------------------------
286	c Monotonic PPM
287	c----------------------------------------------------------------------
288	else if (iadv(iq).eq.16) then
289	c Ss-prg PPM3d de Lin
290	call ppm3d(1,qppm(1,1,iq),
291	s psppm,psppm,
292	s unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
293	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
294	s fill,dum,220.)
295	c---------------------------------------------------------------------
296
297	c---------------------------------------------------------------------
298	c Semi Monotonic PPM
299	c---------------------------------------------------------------------
300	else if (iadv(iq).eq.17) then
301	c Ss-prg PPM3d de Lin
302	call ppm3d(1,qppm(1,1,iq),
303	s psppm,psppm,
304	s unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
305	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
306	s fill,dum,220.)
307	c---------------------------------------------------------------------
308
309	c---------------------------------------------------------------------
310	c Positive Definite PPM
311	c---------------------------------------------------------------------
312	else if (iadv(iq).eq.18) then
313	c Ss-prg PPM3d de Lin
314	call ppm3d(1,qppm(1,1,iq),
315	s psppm,psppm,
316	s unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
317	s iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
318	s fill,dum,220.)
319	c---------------------------------------------------------------------
320	endif
321	enddo
322	c-----------------------------------------------------------------
323	c Ss-prg interface PPM3d-LMDZ.4
324	c-----------------------------------------------------------------
325	call interpost(q(1,1,iq),qppm(1,1,iq))
326	endif
327	c----------------------------------------------------------------------
328
329	c-----------------------------------------------------------------
330	c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
331	c et Nord j=1
332	c-----------------------------------------------------------------
333
334	c call traceurpole(q(1,1,iq),massem)
335
336	c calcul du temps cpu pour un schema donne
337
338
339	end DO
340
341	1234 CONTINUE
342	c$OMP BARRIER
343
344	if (planet_type=="earth") then
345
346	ijb=ij_begin
347	ije=ij_end
348
349	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
350	DO l = 1, llm
351	DO ij = ijb, ije
352	finmasse(ij,l) = p(ij,l) - p(ij,l+1)
353	ENDDO
354	ENDDO
355	c$OMP END DO
356
357	CALL qminimum_loc( q, 2, finmasse )
358
359	endif ! of if (planet_type=="earth")
360
361	RETURN
362	END
363

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

Download in other formats:

Original Format