Context Navigation

vlspltgen_p.F @ 853

Last change on this file since 853 was 774, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 9.6 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	c$OMP MASTER
90	call SetTag(MyRequest1,100)
91	call SetTag(MyRequest2,101)
92	c$OMP END MASTER
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	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
186	& ij_begin,ij_begin+2*iip1-1)
187	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
188	& ij_end-2*iip1+1,ij_end)
189
190	c$OMP MASTER
191	call VTb(VTHallo)
192	call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
193	call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
194	call VTe(VTHallo)
195	c$OMP END MASTER
196	else if (iadv(iq)==14) then
197
198	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
199	& ij_begin,ij_begin+2*iip1-1)
200	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
201	& ij_end-2*iip1+1,ij_end)
202	c$OMP MASTER
203	call VTb(VTHallo)
204	call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
205	call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
206	call VTe(VTHallo)
207	c$OMP END MASTER
208	else
209
210	stop 'vlspltgen_p : schema non parallelise'
211
212	endif
213
214	enddo
215
216
217	c$OMP BARRIER
218	c$OMP MASTER
219	call VTb(VTHallo)
220	call SendRequest(MyRequest1)
221	call VTe(VTHallo)
222	c$OMP END MASTER
223	c$OMP BARRIER
224	do iq=1,nqmx
225
226	if(iadv(iq) == 0) then
227
228	cycle
229
230	else if (iadv(iq)==10) then
231
232	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
233	& ij_begin+2iip1,ij_end-2iip1)
234
235	else if (iadv(iq)==14) then
236
237	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
238	& ij_begin+2iip1,ij_end-2iip1)
239
240	else
241
242	stop 'vlspltgen_p : schema non parallelise'
243
244	endif
245
246	enddo
247	c$OMP BARRIER
248	c$OMP MASTER
249	call VTb(VTHallo)
250	call WaitRecvRequest(MyRequest1)
251	call WaitSendRequest(MyRequest1)
252	call VTe(VTHallo)
253	c$OMP END MASTER
254	c$OMP BARRIER
255
256	do iq=1,nqmx
257
258	if(iadv(iq) == 0) then
259
260	cycle
261
262	else if (iadv(iq)==10) then
263
264	call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
265
266	else if (iadv(iq)==14) then
267
268	call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
269
270	else
271
272	stop 'vlspltgen_p : schema non parallelise'
273
274	endif
275
276	enddo
277
278	do iq=1,nqmx
279
280	if(iadv(iq) == 0) then
281
282	cycle
283
284	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
285
286	c$OMP BARRIER
287	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
288	& ij_begin,ij_begin+2*iip1-1)
289	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
290	& ij_end-2*iip1+1,ij_end)
291	c$OMP BARRIER
292
293	c$OMP MASTER
294	call VTb(VTHallo)
295	call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
296	call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
297	call VTe(VTHallo)
298	c$OMP END MASTER
299	c$OMP BARRIER
300	else
301
302	stop 'vlspltgen_p : schema non parallelise'
303
304	endif
305
306	enddo
307	c$OMP BARRIER
308	c$OMP MASTER
309	call VTb(VTHallo)
310	call SendRequest(MyRequest2)
311	call VTe(VTHallo)
312	c$OMP END MASTER
313	c$OMP BARRIER
314	do iq=1,nqmx
315
316	if(iadv(iq) == 0) then
317
318	cycle
319
320	else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
321	c$OMP BARRIER
322	call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
323	& ij_begin+2iip1,ij_end-2iip1)
324	c$OMP BARRIER
325	else
326
327	stop 'vlspltgen_p : schema non parallelise'
328
329	endif
330
331	enddo
332
333	c$OMP BARRIER
334	c$OMP MASTER
335	call VTb(VTHallo)
336	call WaitRecvRequest(MyRequest2)
337	call WaitSendRequest(MyRequest2)
338	call VTe(VTHallo)
339	c$OMP END MASTER
340	c$OMP BARRIER
341
342
343	do iq=1,nqmx
344
345	if(iadv(iq) == 0) then
346
347	cycle
348
349	else if (iadv(iq)==10) then
350
351	call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
352
353	else if (iadv(iq)==14) then
354
355	call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
356
357	else
358
359	stop 'vlspltgen_p : schema non parallelise'
360
361	endif
362
363	enddo
364
365	do iq=1,nqmx
366
367	if(iadv(iq) == 0) then
368
369	cycle
370
371	else if (iadv(iq)==10) then
372
373	call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
374	& ij_begin,ij_end)
375
376	else if (iadv(iq)==14) then
377
378	call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
379	& ij_begin,ij_end)
380
381	else
382
383	stop 'vlspltgen_p : schema non parallelise'
384
385	endif
386
387	enddo
388
389
390	ijb=ij_begin
391	ije=ij_end
392	c$OMP BARRIER
393
394	DO iq=1,nqmx
395
396	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
397	DO l=1,llm
398	DO ij=ijb,ije
399	c print *,'zq-->',ij,l,iq,zq(ij,l,iq)
400	c print *,'q-->',ij,l,iq,q(ij,l,iq)
401	q(ij,l,iq)=zq(ij,l,iq)
402	ENDDO
403	ENDDO
404	c$OMP END DO NOWAIT
405
406	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
407	DO l=1,llm
408	DO ij=ijb,ije-iip1+1,iip1
409	q(ij+iim,l,iq)=q(ij,l,iq)
410	ENDDO
411	ENDDO
412	c$OMP END DO NOWAIT
413
414	ENDDO
415
416
417	c$OMP BARRIER
418
419	cc$OMP MASTER
420	c call WaitSendRequest(MyRequest1)
421	c call WaitSendRequest(MyRequest2)
422	cc$OMP END MASTER
423	cc$OMP BARRIER
424
425	RETURN
426	END

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

Download in other formats:

Original Format