Context Navigation

vlspltgen_p.F @ 1012

Last change on this file since 1012 was 985, checked in by Laurent Fairhead, 16 years ago
Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.4 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE vlspltgen_p( q,iadv,pente_max,masse,w,pbaru,pbarv,pdt,
5	, p,pk,teta )
6	c
7	c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
8	c
9	c ********************************************************************
10	c Shema d'advection " pseudo amont " .
11	c + test sur humidite specifique: Q advecte< Qsat aval
12	c (F. Codron, 10/99)
13	c ********************************************************************
14	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
15	c
16	c pente_max facteur de limitation des pentes: 2 en general
17	c 0 pour un schema amont
18	c pbaru,pbarv,w flux de masse en u ,v ,w
19	c pdt pas de temps
20	c
21	c teta temperature potentielle, p pression aux interfaces,
22	c pk exner au milieu des couches necessaire pour calculer Qsat
23	c --------------------------------------------------------------------
24	USE parallel
25	USE mod_hallo
26	USE Write_Field_p
27	USE VAMPIR
28	IMPLICIT NONE
29
30	c
31	#include "dimensions.h"
32	#include "paramet.h"
33	#include "logic.h"
34	#include "comvert.h"
35	#include "comconst.h"
36
37	c
38	c Arguments:
39	c ----------
40	INTEGER iadv(nqmx)
41	REAL masse(ip1jmp1,llm),pente_max
42	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
43	REAL q(ip1jmp1,llm,nqmx)
44	REAL w(ip1jmp1,llm),pdt
45	REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
46	c
47	c Local
48	c ---------
49	c
50	INTEGER ij,l
51	c
52	REAL,SAVE :: qsat(ip1jmp1,llm)
53	REAL,SAVE :: zm(ip1jmp1,llm,nqmx)
54	REAL,SAVE :: mu(ip1jmp1,llm)
55	REAL,SAVE :: mv(ip1jm,llm)
56	REAL,SAVE :: mw(ip1jmp1,llm+1)
57	REAL,SAVE :: zq(ip1jmp1,llm,nqmx)
58	REAL zzpbar, zzw
59
60	REAL qmin,qmax
61	DATA qmin,qmax/0.,1.e33/
62
63	c--pour rapport de melange saturant--
64
65	REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
66	REAL ptarg,pdelarg,foeew,zdelta
67	REAL tempe(ip1jmp1)
68	INTEGER ijb,ije,iq
69	type(request) :: MyRequest1
70	type(request) :: MyRequest2
71
72	c fonction psat(T)
73
74	FOEEW ( PTARG,PDELARG ) = EXP (
75	* (R3LES(1.-PDELARG)+R3IESPDELARG) * (PTARG-RTT)
76	* / (PTARG-(R4LES(1.-PDELARG)+R4IESPDELARG)) )
77
78	r2es = 380.11733
79	r3les = 17.269
80	r3ies = 21.875
81	r4les = 35.86
82	r4ies = 7.66
83	retv = 0.6077667
84	rtt = 273.16
85
86	c-- Calcul de Qsat en chaque point
87	c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
88	c pour eviter une exponentielle.
89
90	call SetTag(MyRequest1,100)
91	call SetTag(MyRequest2,101)
92
93
94	ijb=ij_begin-iip1
95	ije=ij_end+iip1
96	if (pole_nord) ijb=ij_begin
97	if (pole_sud) ije=ij_end
98
99	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
100	DO l = 1, llm
101	DO ij = ijb, ije
102	tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
103	ENDDO
104	DO ij = ijb, ije
105	zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
106	play = 0.5*(p(ij,l)+p(ij,l+1))
107	qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
108	qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
109	ENDDO
110	ENDDO
111	c$OMP END DO NOWAIT
112	c PRINT*,'Debut vlsplt version debug sans vlyqs'
113
114	zzpbar = 0.5 * pdt
115	zzw = pdt
116
117	ijb=ij_begin
118	ije=ij_end
119	if (pole_nord) ijb=ijb+iip1
120	if (pole_sud) ije=ije-iip1
121
122	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
123	DO l=1,llm
124	DO ij = ijb,ije
125	mu(ij,l)=pbaru(ij,l) * zzpbar
126	ENDDO
127	ENDDO
128	c$OMP END DO NOWAIT
129
130	ijb=ij_begin-iip1
131	ije=ij_end
132	if (pole_nord) ijb=ij_begin
133	if (pole_sud) ije=ij_end-iip1
134
135	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
136	DO l=1,llm
137	DO ij=ijb,ije
138	mv(ij,l)=pbarv(ij,l) * zzpbar
139	ENDDO
140	ENDDO
141	c$OMP END DO NOWAIT
142
143	ijb=ij_begin
144	ije=ij_end
145
146	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
147	DO l=1,llm
148	DO ij=ijb,ije
149	mw(ij,l)=w(ij,l) * zzw
150	ENDDO
151	ENDDO
152	c$OMP END DO NOWAIT
153
154	c$OMP MASTER
155	DO ij=ijb,ije
156	mw(ij,llm+1)=0.
157	ENDDO
158	c$OMP END MASTER
159
160	c CALL SCOPY(ijp1llm,q,1,zq,1)
161	c CALL SCOPY(ijp1llm,masse,1,zm,1)
162
163	ijb=ij_begin
164	ije=ij_end
165
166	DO iq=1,nqmx
167	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168	DO l=1,llm
169	zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
170	zm(ijb:ije,l,iq)=masse(ijb:ije,l)
171	ENDDO
172	c$OMP END DO NOWAIT
173	ENDDO
174
175
176	c$OMP BARRIER
177	DO iq=1,nqmx
178
179	if(iadv(iq) == 0) then
180
181	cycle
182
183	else if (iadv(iq)==10) then
184
185	#ifdef _ADV_HALO
186	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
187	& ij_begin,ij_begin+2*iip1-1)
188	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
189	& ij_end-2*iip1+1,ij_end)
190	#else
191	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
192	& ij_begin,ij_end)
193	#endif
194
195	c$OMP MASTER
196	call VTb(VTHallo)
197	c$OMP END MASTER
198	call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
199	call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
200
201	c$OMP MASTER
202	call VTe(VTHallo)
203	c$OMP END MASTER
204	else if (iadv(iq)==14) then
205
206	#ifdef _ADV_HALO
207	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
208	& ij_begin,ij_begin+2*iip1-1)
209	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
210	& ij_end-2*iip1+1,ij_end)
211	#else
212
213	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
214	& ij_begin,ij_end)
215	#endif
216
217	c$OMP MASTER
218	call VTb(VTHallo)
219	c$OMP END MASTER
220
221	call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
222	call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
223
224	c$OMP MASTER
225	call VTe(VTHallo)
226	c$OMP END MASTER
227	else
228
229	stop 'vlspltgen_p : schema non parallelise'
230
231	endif
232
233	enddo
234
235
236	c$OMP BARRIER
237	c$OMP MASTER
238	call VTb(VTHallo)
239	c$OMP END MASTER
240
241	call SendRequest(MyRequest1)
242
243	c$OMP MASTER
244	call VTe(VTHallo)
245	c$OMP END MASTER
246	c$OMP BARRIER
247	do iq=1,nqmx
248
249	if(iadv(iq) == 0) then
250
251	cycle
252
253	else if (iadv(iq)==10) then
254
255	#ifdef _ADV_HALLO
256	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
257	& ij_begin+2iip1,ij_end-2iip1)
258	#endif
259	else if (iadv(iq)==14) then
260	#ifdef _ADV_HALLO
261	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
262	& ij_begin+2iip1,ij_end-2iip1)
263	#endif
264	else
265
266	stop 'vlspltgen_p : schema non parallelise'
267
268	endif
269
270	enddo
271	c$OMP BARRIER
272	c$OMP MASTER
273	call VTb(VTHallo)
274	c$OMP END MASTER
275
276	! call WaitRecvRequest(MyRequest1)
277	! call WaitSendRequest(MyRequest1)
278	c$OMP BARRIER
279	call WaitRequest(MyRequest1)
280
281
282	c$OMP MASTER
283	call VTe(VTHallo)
284	c$OMP END MASTER
285	c$OMP BARRIER
286
287	do iq=1,nqmx
288
289	if(iadv(iq) == 0) then
290
291	cycle
292
293	else if (iadv(iq)==10) then
294
295	call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
296
297	else if (iadv(iq)==14) then
298
299	call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
300
301	else
302
303	stop 'vlspltgen_p : schema non parallelise'
304
305	endif
306
307	enddo
308
309
310	do iq=1,nqmx
311
312	if(iadv(iq) == 0) then
313
314	cycle
315
316	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
317
318	c$OMP BARRIER
319	#ifdef _ADV_HALLO
320	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
321	& ij_begin,ij_begin+2*iip1-1)
322	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
323	& ij_end-2*iip1+1,ij_end)
324	#else
325	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
326	& ij_begin,ij_end)
327	#endif
328	c$OMP BARRIER
329
330	c$OMP MASTER
331	call VTb(VTHallo)
332	c$OMP END MASTER
333
334	call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
335	call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
336
337	c$OMP MASTER
338	call VTe(VTHallo)
339	c$OMP END MASTER
340	c$OMP BARRIER
341	else
342
343	stop 'vlspltgen_p : schema non parallelise'
344
345	endif
346
347	enddo
348	c$OMP BARRIER
349
350	c$OMP MASTER
351	call VTb(VTHallo)
352	c$OMP END MASTER
353
354	call SendRequest(MyRequest2)
355
356	c$OMP MASTER
357	call VTe(VTHallo)
358	c$OMP END MASTER
359
360	c$OMP BARRIER
361	do iq=1,nqmx
362
363	if(iadv(iq) == 0) then
364
365	cycle
366
367	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
368	c$OMP BARRIER
369
370	#ifdef _ADV_HALLO
371	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
372	& ij_begin+2iip1,ij_end-2iip1)
373	#endif
374
375	c$OMP BARRIER
376	else
377
378	stop 'vlspltgen_p : schema non parallelise'
379
380	endif
381
382	enddo
383
384	c$OMP BARRIER
385	c$OMP MASTER
386	call VTb(VTHallo)
387	c$OMP END MASTER
388
389	! call WaitRecvRequest(MyRequest2)
390	! call WaitSendRequest(MyRequest2)
391	c$OMP BARRIER
392	CALL WaitRequest(MyRequest2)
393
394	c$OMP MASTER
395	call VTe(VTHallo)
396	c$OMP END MASTER
397	c$OMP BARRIER
398
399
400	do iq=1,nqmx
401
402	if(iadv(iq) == 0) then
403
404	cycle
405
406	else if (iadv(iq)==10) then
407
408	call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
409
410	else if (iadv(iq)==14) then
411
412	call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
413
414	else
415
416	stop 'vlspltgen_p : schema non parallelise'
417
418	endif
419
420	enddo
421
422	do iq=1,nqmx
423
424	if(iadv(iq) == 0) then
425
426	cycle
427
428	else if (iadv(iq)==10) then
429
430	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
431	& ij_begin,ij_end)
432
433	else if (iadv(iq)==14) then
434
435	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
436	& ij_begin,ij_end)
437
438	else
439
440	stop 'vlspltgen_p : schema non parallelise'
441
442	endif
443
444	enddo
445
446
447	ijb=ij_begin
448	ije=ij_end
449	c$OMP BARRIER
450
451
452	DO iq=1,nqmx
453
454	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
455	DO l=1,llm
456	DO ij=ijb,ije
457	c print *,'zq-->',ij,l,iq,zq(ij,l,iq)
458	c print *,'q-->',ij,l,iq,q(ij,l,iq)
459	q(ij,l,iq)=zq(ij,l,iq)
460	ENDDO
461	ENDDO
462	c$OMP END DO NOWAIT
463
464	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
465	DO l=1,llm
466	DO ij=ijb,ije-iip1+1,iip1
467	q(ij+iim,l,iq)=q(ij,l,iq)
468	ENDDO
469	ENDDO
470	c$OMP END DO NOWAIT
471
472	ENDDO
473
474
475	c$OMP BARRIER
476
477	cc$OMP MASTER
478	c call WaitSendRequest(MyRequest1)
479	c call WaitSendRequest(MyRequest2)
480	cc$OMP END MASTER
481	cc$OMP BARRIER
482
483	RETURN
484	END

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

Download in other formats:

Original Format