Context Navigation

vlsplt_p.F @ 1114

Last change on this file since 1114 was 764, checked in by Laurent Fairhead, 17 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.8 KB

Line
1	SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
2	c
3	c Auteurs: P.Le Van, F.Hourdin, F.Forget
4	c
5	c ********************************************************************
6	c Shema d'advection " pseudo amont " .
7	c ********************************************************************
8	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
9	c
10	c pente_max facteur de limitation des pentes: 2 en general
11	c 0 pour un schema amont
12	c pbaru,pbarv,w flux de masse en u ,v ,w
13	c pdt pas de temps
14	c
15	c --------------------------------------------------------------------
16	USE parallel
17	USE mod_hallo
18	USE Vampir
19	IMPLICIT NONE
20	c
21	#include "dimensions.h"
22	#include "paramet.h"
23	#include "logic.h"
24	#include "comvert.h"
25	#include "comconst.h"
26
27	c
28	c Arguments:
29	c ----------
30	REAL masse(ip1jmp1,llm),pente_max
31	c REAL masse(iip1,jjp1,llm),pente_max
32	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
33	REAL q(ip1jmp1,llm)
34	c REAL q(iip1,jjp1,llm)
35	REAL w(ip1jmp1,llm),pdt
36	c
37	c Local
38	c ---------
39	c
40	INTEGER i,ij,l,j,ii
41	INTEGER ijlqmin,iqmin,jqmin,lqmin
42	c
43	REAL zm(ip1jmp1,llm),newmasse
44	REAL mu(ip1jmp1,llm)
45	REAL mv(ip1jm,llm)
46	REAL mw(ip1jmp1,llm+1)
47	REAL zq(ip1jmp1,llm),zz
48	REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
49	REAL second,temps0,temps1,temps2,temps3
50	REAL ztemps1,ztemps2,ztemps3
51	REAL zzpbar, zzw
52	LOGICAL testcpu
53	SAVE testcpu
54	SAVE temps1,temps2,temps3
55	INTEGER iminn,imaxx
56
57	REAL qmin,qmax
58	DATA qmin,qmax/0.,1.e33/
59	DATA testcpu/.false./
60	DATA temps1,temps2,temps3/0.,0.,0./
61	INTEGER ijb,ije
62	type(request) :: MyRequest1
63	type(request) :: MyRequest2
64
65	call SetTag(MyRequest1,100)
66	call SetTag(MyRequest2,101)
67
68	zzpbar = 0.5 * pdt
69	zzw = pdt
70
71	ijb=ij_begin
72	ije=ij_end
73	if (pole_nord) ijb=ijb+iip1
74	if (pole_sud) ije=ije-iip1
75
76	DO l=1,llm
77	DO ij = ijb,ije
78	mu(ij,l)=pbaru(ij,l) * zzpbar
79	ENDDO
80	ENDDO
81
82	ijb=ij_begin-iip1
83	ije=ij_end
84	if (pole_nord) ijb=ij_begin
85	if (pole_sud) ije=ij_end-iip1
86
87	DO l=1,llm
88	DO ij=ijb,ije
89	mv(ij,l)=pbarv(ij,l) * zzpbar
90	ENDDO
91	ENDDO
92
93	ijb=ij_begin
94	ije=ij_end
95
96	DO l=1,llm
97	DO ij=ijb,ije
98	mw(ij,l)=w(ij,l) * zzw
99	ENDDO
100	ENDDO
101
102	DO ij=ijb,ije
103	mw(ij,llm+1)=0.
104	ENDDO
105
106	c CALL SCOPY(ijp1llm,q,1,zq,1)
107	c CALL SCOPY(ijp1llm,masse,1,zm,1)
108
109	ijb=ij_begin
110	ije=ij_end
111	zq(ijb:ije,:)=q(ijb:ije,:)
112	zm(ijb:ije,:)=masse(ijb:ije,:)
113
114
115	c print*,'Entree vlx1'
116	c call minmaxq(zq,qmin,qmax,'avant vlx ')
117	call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1)
118	call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end)
119	call VTb(VTHallo)
120	call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
121	call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest1)
122	call SendRequest(MyRequest1)
123	call VTe(VTHallo)
124	call vlx_p(zq,pente_max,zm,mu,ij_begin+2iip1,ij_end-2iip1)
125	c call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
126	call VTb(VTHallo)
127	call WaitRecvRequest(MyRequest1)
128	call VTe(VTHallo)
129
130
131	c print*,'Sortie vlx1'
132	c call minmaxq(zq,qmin,qmax,'apres vlx1 ')
133
134	c print*,'Entree vly1'
135	c call exchange_hallo(zq,ip1jmp1,llm,2,2)
136	c call exchange_hallo(zm,ip1jmp1,llm,1,1)
137
138	call vly_p(zq,pente_max,zm,mv)
139	c call minmaxq(zq,qmin,qmax,'apres vly1 ')
140	c print*,'Sortie vly1'
141	call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
142	call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
143	call VTb(VTHallo)
144	call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
145	call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest2)
146	call SendRequest(MyRequest2)
147	call VTe(VTHallo)
148	call vlz_p(zq,pente_max,zm,mw,ij_begin+2iip1,ij_end-2iip1)
149	call VTb(VTHallo)
150	call WaitRecvRequest(MyRequest2)
151
152	call VTe(VTHallo)
153
154	c call minmaxq(zq,qmin,qmax,'apres vlz ')
155
156
157
158
159	call vly_p(zq,pente_max,zm,mv)
160	c call minmaxq(zq,qmin,qmax,'apres vly ')
161
162
163	call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
164	c call minmaxq(zq,qmin,qmax,'apres vlx2 ')
165
166
167	ijb=ij_begin
168	ije=ij_end
169
170	DO l=1,llm
171	DO ij=ijb,ije
172	q(ij,l)=zq(ij,l)
173	ENDDO
174	ENDDO
175
176
177	DO l=1,llm
178	DO ij=ijb,ije-iip1+1,iip1
179	q(ij+iim,l)=q(ij,l)
180	ENDDO
181	ENDDO
182
183	call WaitSendRequest(MyRequest1)
184	call WaitSendRequest(MyRequest2)
185
186	RETURN
187	END
188
189
190	SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x)
191
192	c Auteurs: P.Le Van, F.Hourdin, F.Forget
193	c
194	c ********************************************************************
195	c Shema d'advection " pseudo amont " .
196	c ********************************************************************
197	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
198	c
199	c
200	c --------------------------------------------------------------------
201	USE Parallel
202	IMPLICIT NONE
203	c
204	#include "dimensions.h"
205	#include "paramet.h"
206	#include "logic.h"
207	#include "comvert.h"
208	#include "comconst.h"
209	c
210	c
211	c Arguments:
212	c ----------
213	REAL masse(ip1jmp1,llm),pente_max
214	REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
215	REAL q(ip1jmp1,llm)
216	REAL w(ip1jmp1,llm)
217	c
218	c Local
219	c ---------
220	c
221	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
222	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
223	c
224	REAL new_m,zu_m,zdum(ip1jmp1,llm)
225	REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
226	REAL zz(ip1jmp1)
227	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
228	REAL u_mq(ip1jmp1,llm)
229
230	Logical extremum
231
232	REAL SSUM
233	EXTERNAL SSUM
234
235	REAL z1,z2,z3
236
237	INTEGER ijb,ije,ijb_x,ije_x
238
239	c calcul de la pente a droite et a gauche de la maille
240
241	ijb=ijb_x
242	ije=ije_x
243
244	if (pole_nord.and.ijb==1) ijb=ijb+iip1
245	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
246
247	IF (pente_max.gt.-1.e-5) THEN
248	c IF (pente_max.gt.10) THEN
249
250	c calcul des pentes avec limitation, Van Leer scheme I:
251	c -----------------------------------------------------
252
253	c calcul de la pente aux points u
254	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
255	DO l = 1, llm
256
257	DO ij=ijb,ije-1
258	dxqu(ij)=q(ij+1,l)-q(ij,l)
259	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
260	c sigu(ij)=u_m(ij,l)/masse(ij,l)
261	ENDDO
262	DO ij=ijb+iip1-1,ije,iip1
263	dxqu(ij)=dxqu(ij-iim)
264	c sigu(ij)=sigu(ij-iim)
265	ENDDO
266
267	DO ij=ijb,ije
268	adxqu(ij)=abs(dxqu(ij))
269	ENDDO
270
271	c calcul de la pente maximum dans la maille en valeur absolue
272
273	DO ij=ijb+1,ije
274	dxqmax(ij,l)=pente_max*
275	, min(adxqu(ij-1),adxqu(ij))
276	c limitation subtile
277	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
278
279
280	ENDDO
281
282	DO ij=ijb+iip1-1,ije,iip1
283	dxqmax(ij-iim,l)=dxqmax(ij,l)
284	ENDDO
285
286	DO ij=ijb+1,ije
287	#ifdef CRAY
288	dxq(ij,l)=
289	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
290	#else
291	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
292	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
293	ELSE
294	c extremum local
295	dxq(ij,l)=0.
296	ENDIF
297	#endif
298	dxq(ij,l)=0.5*dxq(ij,l)
299	dxq(ij,l)=
300	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
301	ENDDO
302
303	ENDDO ! l=1,llm
304	c$OMP END DO NOWAIT
305	c print*,'Ok calcul des pentes'
306
307	ELSE ! (pente_max.lt.-1.e-5)
308
309	c Pentes produits:
310	c ----------------
311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312	DO l = 1, llm
313	DO ij=ijb,ije-1
314	dxqu(ij)=q(ij+1,l)-q(ij,l)
315	ENDDO
316	DO ij=ijb+iip1-1,ije,iip1
317	dxqu(ij)=dxqu(ij-iim)
318	ENDDO
319
320	DO ij=ijb+1,ije
321	zz(ij)=dxqu(ij-1)*dxqu(ij)
322	zz(ij)=zz(ij)+zz(ij)
323	IF(zz(ij).gt.0) THEN
324	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
325	ELSE
326	c extremum local
327	dxq(ij,l)=0.
328	ENDIF
329	ENDDO
330
331	ENDDO
332	c$OMP END DO NOWAIT
333	ENDIF ! (pente_max.lt.-1.e-5)
334
335	c bouclage de la pente en iip1:
336	c -----------------------------
337	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
338	DO l=1,llm
339	DO ij=ijb+iip1-1,ije,iip1
340	dxq(ij-iim,l)=dxq(ij,l)
341	ENDDO
342	DO ij=ijb,ije
343	iadvplus(ij,l)=0
344	ENDDO
345
346	ENDDO
347	c$OMP END DO NOWAIT
348	c print*,'Bouclage en iip1'
349
350	c calcul des flux a gauche et a droite
351
352	#ifdef CRAY
353	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
354	DO l=1,llm
355	DO ij=ijb,ije-1
356	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
357	, 1.+u_m(ij,l)/masse(ij+1,l),
358	, u_m(ij,l))
359	zdum(ij,l)=0.5*zdum(ij,l)
360	u_mq(ij,l)=cvmgp(
361	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
362	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
363	, u_m(ij,l))
364	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
365	ENDDO
366	ENDDO
367	c$OMP END DO NOWAIT
368	#else
369	c on cumule le flux correspondant a toutes les mailles dont la masse
370	c au travers de la paroi pENDant le pas de temps.
371	c print*,'Cumule ....'
372	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
373	DO l=1,llm
374	DO ij=ijb,ije-1
375	c print*,'masse(',ij,')=',masse(ij,l)
376	IF (u_m(ij,l).gt.0.) THEN
377	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
378	u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
379	ELSE
380	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
381	u_mq(ij,l)=u_m(ij,l)(q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
382	ENDIF
383	ENDDO
384	ENDDO
385	c$OMP END DO NOWAIT
386	#endif
387	c stop
388
389	c go to 9999
390	c detection des points ou on advecte plus que la masse de la
391	c maille
392	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
393	DO l=1,llm
394	DO ij=ijb,ije-1
395	IF(zdum(ij,l).lt.0) THEN
396	iadvplus(ij,l)=1
397	u_mq(ij,l)=0.
398	ENDIF
399	ENDDO
400	ENDDO
401	c$OMP END DO NOWAIT
402	c print*,'Ok test 1'
403	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
404	DO l=1,llm
405	DO ij=ijb+iip1-1,ije,iip1
406	iadvplus(ij,l)=iadvplus(ij-iim,l)
407	ENDDO
408	ENDDO
409	c$OMP END DO NOWAIT
410	c print*,'Ok test 2'
411
412
413	c traitement special pour le cas ou on advecte en longitude plus que le
414	c contenu de la maille.
415	c cette partie est mal vectorisee.
416
417	c calcul du nombre de maille sur lequel on advecte plus que la maille.
418
419	n0=0
420	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
421	DO l=1,llm
422	nl(l)=0
423	DO ij=ijb,ije
424	nl(l)=nl(l)+iadvplus(ij,l)
425	ENDDO
426	n0=n0+nl(l)
427	ENDDO
428	c$OMP END DO NOWAIT
429	cym IF(n0.gt.1) THEN
430	cym IF(n0.gt.0) THEN
431
432	c PRINT*,'Nombre de points pour lesquels on advect plus que le'
433	c & ,'contenu de la maille : ',n0
434	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
435	DO l=1,llm
436	IF(nl(l).gt.0) THEN
437	iju=0
438	c indicage des mailles concernees par le traitement special
439	DO ij=ijb,ije
440	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
441	iju=iju+1
442	indu(iju)=ij
443	ENDIF
444	ENDDO
445	niju=iju
446	c PRINT*,'niju,nl',niju,nl(l)
447
448	c traitement des mailles
449	DO iju=1,niju
450	ij=indu(iju)
451	j=(ij-1)/iip1+1
452	zu_m=u_m(ij,l)
453	u_mq(ij,l)=0.
454	IF(zu_m.gt.0.) THEN
455	ijq=ij
456	i=ijq-(j-1)*iip1
457	c accumulation pour les mailles completements advectees
458	do while(zu_m.gt.masse(ijq,l))
459	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
460	zu_m=zu_m-masse(ijq,l)
461	i=mod(i-2+iim,iim)+1
462	ijq=(j-1)*iip1+i
463	ENDDO
464	c ajout de la maille non completement advectee
465	u_mq(ij,l)=u_mq(ij,l)+zu_m*
466	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
467	ELSE
468	ijq=ij+1
469	i=ijq-(j-1)*iip1
470	c accumulation pour les mailles completements advectees
471	do while(-zu_m.gt.masse(ijq,l))
472	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
473	zu_m=zu_m+masse(ijq,l)
474	i=mod(i,iim)+1
475	ijq=(j-1)*iip1+i
476	ENDDO
477	c ajout de la maille non completement advectee
478	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
479	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
480	ENDIF
481	ENDDO
482	ENDIF
483	ENDDO
484	c$OMP END DO NOWAIT
485	cym ENDIF ! n0.gt.0
486	9999 continue
487
488
489	c bouclage en latitude
490	c print*,'Avant bouclage en latitude'
491	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
492	DO l=1,llm
493	DO ij=ijb+iip1-1,ije,iip1
494	u_mq(ij,l)=u_mq(ij-iim,l)
495	ENDDO
496	ENDDO
497	c$OMP END DO NOWAIT
498
499	c calcul des tENDances
500	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
501	DO l=1,llm
502	DO ij=ijb+1,ije
503	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
504	q(ij,l)=(q(ij,l)*masse(ij,l)+
505	& u_mq(ij-1,l)-u_mq(ij,l))
506	& /new_m
507	masse(ij,l)=new_m
508	ENDDO
509	c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
510	DO ij=ijb+iip1-1,ije,iip1
511	q(ij-iim,l)=q(ij,l)
512	masse(ij-iim,l)=masse(ij,l)
513	ENDDO
514	ENDDO
515	c$OMP END DO NOWAIT
516	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
517	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
518
519
520	RETURN
521	END
522
523
524	SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
525	c
526	c Auteurs: P.Le Van, F.Hourdin, F.Forget
527	c
528	c ********************************************************************
529	c Shema d'advection " pseudo amont " .
530	c ********************************************************************
531	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
532	c dq sont des arguments de sortie pour le s-pg ....
533	c
534	c
535	c --------------------------------------------------------------------
536	USE parallel
537	IMPLICIT NONE
538	c
539	#include "dimensions.h"
540	#include "paramet.h"
541	#include "logic.h"
542	#include "comvert.h"
543	#include "comconst.h"
544	#include "comgeom.h"
545	c
546	c
547	c Arguments:
548	c ----------
549	REAL masse(ip1jmp1,llm),pente_max
550	REAL masse_adv_v( ip1jm,llm)
551	REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
552	c
553	c Local
554	c ---------
555	c
556	INTEGER i,ij,l
557	c
558	REAL airej2,airejjm,airescb(iim),airesch(iim)
559	REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
560	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
561	REAL qbyv(ip1jm,llm)
562
563	REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
564	c REAL newq,oldmasse
565	Logical extremum,first,testcpu
566	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
567	SAVE temps0,temps1,temps2,temps3,temps4,temps5
568	c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
569	SAVE first,testcpu
570	c$OMP THREADPRIVATE(first,testcpu)
571
572	REAL convpn,convps,convmpn,convmps
573	real massepn,masseps,qpn,qps
574	REAL sinlon(iip1),sinlondlon(iip1)
575	REAL coslon(iip1),coslondlon(iip1)
576	SAVE sinlon,coslon,sinlondlon,coslondlon
577	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
578	SAVE airej2,airejjm
579	c$OMP THREADPRIVATE(airej2,airejjm)
580	c
581	c
582	REAL SSUM
583	EXTERNAL SSUM
584
585	DATA first,testcpu/.true.,.false./
586	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
587	INTEGER ijb,ije
588
589	IF(first) THEN
590	c PRINT*,'Shema Amont nouveau appele dans Vanleer '
591	first=.false.
592	do i=2,iip1
593	coslon(i)=cos(rlonv(i))
594	sinlon(i)=sin(rlonv(i))
595	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
596	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
597	ENDDO
598	coslon(1)=coslon(iip1)
599	coslondlon(1)=coslondlon(iip1)
600	sinlon(1)=sinlon(iip1)
601	sinlondlon(1)=sinlondlon(iip1)
602	airej2 = SSUM( iim, aire(iip2), 1 )
603	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
604	ENDIF
605
606	c
607	c PRINT*,'CALCUL EN LATITUDE'
608
609	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
610	DO l = 1, llm
611	c
612	c --------------------------------
613	c CALCUL EN LATITUDE
614	c --------------------------------
615
616	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
617	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
618	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
619
620	if (pole_nord) then
621	DO i = 1, iim
622	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
623	ENDDO
624	qpns = SSUM( iim, airescb ,1 ) / airej2
625	endif
626
627	if (pole_sud) then
628	DO i = 1, iim
629	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
630	ENDDO
631	qpsn = SSUM( iim, airesch ,1 ) / airejjm
632	endif
633
634
635
636	c calcul des pentes aux points v
637
638	ijb=ij_begin-2*iip1
639	ije=ij_end+iip1
640	if (pole_nord) ijb=ij_begin
641	if (pole_sud) ije=ij_end-iip1
642
643	DO ij=ijb,ije
644	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
645	adyqv(ij)=abs(dyqv(ij))
646	ENDDO
647
648	c calcul des pentes aux points scalaires
649	ijb=ij_begin-iip1
650	ije=ij_end+iip1
651	if (pole_nord) ijb=ij_begin+iip1
652	if (pole_sud) ije=ij_end-iip1
653
654	DO ij=ijb,ije
655	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
656	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
657	dyqmax(ij)=pente_max*dyqmax(ij)
658	ENDDO
659
660	c calcul des pentes aux poles
661	IF (pole_nord) THEN
662	DO ij=1,iip1
663	dyq(ij,l)=qpns-q(ij+iip1,l)
664	ENDDO
665
666	dyn1=0.
667	dyn2=0.
668	DO ij=1,iim
669	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
670	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
671	ENDDO
672	DO ij=1,iip1
673	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
674	ENDDO
675
676	DO ij=1,iip1
677	dyq(ij,l)=0.
678	ENDDO
679	c ym tout cela ne sert pas a grand chose
680	ENDIF
681
682	IF (pole_sud) THEN
683
684	DO ij=1,iip1
685	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
686	ENDDO
687
688	dys1=0.
689	dys2=0.
690
691	DO ij=1,iim
692	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
693	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
694	ENDDO
695
696	DO ij=1,iip1
697	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
698	ENDDO
699
700	DO ij=1,iip1
701	dyq(ip1jm+ij,l)=0.
702	ENDDO
703	c ym tout cela ne sert pas a grand chose
704	ENDIF
705
706	c filtrage de la derivee
707
708	c calcul des pentes limites aux poles
709	c ym partie inutile
710	c goto 8888
711	c fn=1.
712	c fs=1.
713	c DO ij=1,iim
714	c IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
715	c fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
716	c ENDIF
717	c IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
718	c fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
719	c ENDIF
720	c ENDDO
721	c DO ij=1,iip1
722	c dyq(ij,l)=fn*dyq(ij,l)
723	c dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
724	c ENDDO
725	c 8888 continue
726
727
728	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
729	C En memoire de dIFferents tests sur la
730	C limitation des pentes aux poles.
731	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
732	C PRINT*,dyq(1)
733	C PRINT*,dyqv(iip1+1)
734	C apn=abs(dyq(1)/dyqv(iip1+1))
735	C PRINT*,dyq(ip1jm+1)
736	C PRINT*,dyqv(ip1jm-iip1+1)
737	C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
738	C DO ij=2,iim
739	C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
740	C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
741	C ENDDO
742	C apn=min(pente_max/apn,1.)
743	C aps=min(pente_max/aps,1.)
744	C
745	C
746	C cas ou on a un extremum au pole
747	C
748	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
749	C & apn=0.
750	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
751	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
752	C & aps=0.
753	C
754	C limitation des pentes aux poles
755	C DO ij=1,iip1
756	C dyq(ij)=apn*dyq(ij)
757	C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
758	C ENDDO
759	C
760	C test
761	C DO ij=1,iip1
762	C dyq(iip1+ij)=0.
763	C dyq(ip1jm+ij-iip1)=0.
764	C ENDDO
765	C DO ij=1,ip1jmp1
766	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
767	C ENDDO
768	C
769	C changement 10 07 96
770	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
771	C & THEN
772	C DO ij=1,iip1
773	C dyqmax(ij)=0.
774	C ENDDO
775	C ELSE
776	C DO ij=1,iip1
777	C dyqmax(ij)=pente_max*abs(dyqv(ij))
778	C ENDDO
779	C ENDIF
780	C
781	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
782	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
783	C &THEN
784	C DO ij=ip1jm+1,ip1jmp1
785	C dyqmax(ij)=0.
786	C ENDDO
787	C ELSE
788	C DO ij=ip1jm+1,ip1jmp1
789	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
790	C ENDDO
791	C ENDIF
792	C fin changement 10 07 96
793	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
794
795	c calcul des pentes limitees
796	ijb=ij_begin-iip1
797	ije=ij_end+iip1
798	if (pole_nord) ijb=ij_begin+iip1
799	if (pole_sud) ije=ij_end-iip1
800
801	DO ij=ijb,ije
802	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
803	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
804	ELSE
805	dyq(ij,l)=0.
806	ENDIF
807	ENDDO
808
809	ENDDO
810	c$OMP END DO NOWAIT
811
812	ijb=ij_begin-iip1
813	ije=ij_end
814	if (pole_nord) ijb=ij_begin
815	if (pole_sud) ije=ij_end-iip1
816
817	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
818	DO l=1,llm
819	DO ij=ijb,ije
820	IF(masse_adv_v(ij,l).gt.0) THEN
821	qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
822	, 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
823	ELSE
824	qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
825	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
826	ENDIF
827	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
828	ENDDO
829	ENDDO
830	c$OMP END DO NOWAIT
831
832	ijb=ij_begin
833	ije=ij_end
834	if (pole_nord) ijb=ij_begin+iip1
835	if (pole_sud) ije=ij_end-iip1
836
837	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
838	DO l=1,llm
839	DO ij=ijb,ije
840	newmasse=masse(ij,l)
841	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
842
843	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
844	& /newmasse
845	masse(ij,l)=newmasse
846	ENDDO
847	c.-. ancienne version
848	c convpn=SSUM(iim,qbyv(1,l),1)/apoln
849	c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
850	if (pole_nord) then
851	convpn=SSUM(iim,qbyv(1,l),1)
852	convmpn=ssum(iim,masse_adv_v(1,l),1)
853	massepn=ssum(iim,masse(1,l),1)
854	qpn=0.
855	do ij=1,iim
856	qpn=qpn+masse(ij,l)*q(ij,l)
857	enddo
858	qpn=(qpn+convpn)/(massepn+convmpn)
859	do ij=1,iip1
860	q(ij,l)=qpn
861	enddo
862	endif
863
864	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
865	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
866
867	if (pole_sud) then
868
869	convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
870	convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
871	masseps=ssum(iim, masse(ip1jm+1,l),1)
872	qps=0.
873	do ij = ip1jm+1,ip1jmp1-1
874	qps=qps+masse(ij,l)*q(ij,l)
875	enddo
876	qps=(qps+convps)/(masseps+convmps)
877	do ij=ip1jm+1,ip1jmp1
878	q(ij,l)=qps
879	enddo
880	endif
881	c.-. fin ancienne version
882
883	c._. nouvelle version
884	c convpn=SSUM(iim,qbyv(1,l),1)
885	c convmpn=ssum(iim,masse_adv_v(1,l),1)
886	c oldmasse=ssum(iim,masse(1,l),1)
887	c newmasse=oldmasse+convmpn
888	c newq=(q(1,l)*oldmasse+convpn)/newmasse
889	c newmasse=newmasse/apoln
890	c DO ij = 1,iip1
891	c q(ij,l)=newq
892	c masse(ij,l)=newmasse*aire(ij)
893	c ENDDO
894	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
895	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
896	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
897	c newmasse=oldmasse+convmps
898	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
899	c newmasse=newmasse/apols
900	c DO ij = ip1jm+1,ip1jmp1
901	c q(ij,l)=newq
902	c masse(ij,l)=newmasse*aire(ij)
903	c ENDDO
904	c._. fin nouvelle version
905	ENDDO
906	c$OMP END DO NOWAIT
907
908	RETURN
909	END
910
911
912
913	SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
914	c
915	c Auteurs: P.Le Van, F.Hourdin, F.Forget
916	c
917	c ********************************************************************
918	c Shema d'advection " pseudo amont " .
919	c ********************************************************************
920	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
921	c dq sont des arguments de sortie pour le s-pg ....
922	c
923	c
924	c --------------------------------------------------------------------
925	USE Parallel
926	IMPLICIT NONE
927	c
928	#include "dimensions.h"
929	#include "paramet.h"
930	#include "logic.h"
931	#include "comvert.h"
932	#include "comconst.h"
933	c
934	c
935	c Arguments:
936	c ----------
937	REAL masse(ip1jmp1,llm),pente_max
938	REAL q(ip1jmp1,llm)
939	REAL w(ip1jmp1,llm+1)
940	c
941	c Local
942	c ---------
943	c
944	INTEGER i,ij,l,j,ii
945	c
946	REAL,SAVE :: wq(ip1jmp1,llm+1)
947	REAL newmasse
948
949	REAL,SAVE :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm)
950	REAL dzqmax
951	REAL sigw
952
953	LOGICAL testcpu
954	SAVE testcpu
955	c$OMP THREADPRIVATE(testcpu)
956	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
957	SAVE temps0,temps1,temps2,temps3,temps4,temps5
958	c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
959
960	REAL SSUM
961	EXTERNAL SSUM
962
963	DATA testcpu/.false./
964	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
965	INTEGER ijb,ije,ijb_x,ije_x
966	c On oriente tout dans le sens de la pression c'est a dire dans le
967	c sens de W
968
969	#ifdef BIDON
970	IF(testcpu) THEN
971	temps0=second(0.)
972	ENDIF
973	#endif
974
975	ijb=ijb_x
976	ije=ije_x
977
978	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
979	DO l=2,llm
980	DO ij=ijb,ije
981	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
982	adzqw(ij,l)=abs(dzqw(ij,l))
983	ENDDO
984	ENDDO
985	c$OMP END DO
986
987	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
988	DO l=2,llm-1
989	DO ij=ijb,ije
990	#ifdef CRAY
991	dzq(ij,l)=0.5*
992	, cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
993	#else
994	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
995	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
996	ELSE
997	dzq(ij,l)=0.
998	ENDIF
999	#endif
1000	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
1001	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
1002	ENDDO
1003	ENDDO
1004	c$OMP END DO NOWAIT
1005
1006	c$OMP MASTER
1007	DO ij=ijb,ije
1008	dzq(ij,1)=0.
1009	dzq(ij,llm)=0.
1010	ENDDO
1011	c$OMP END MASTER
1012	c$OMP BARRIER
1013	#ifdef BIDON
1014	IF(testcpu) THEN
1015	temps1=temps1+second(0.)-temps0
1016	ENDIF
1017	#endif
1018	c ---------------------------------------------------------------
1019	c .... calcul des termes d'advection verticale .......
1020	c ---------------------------------------------------------------
1021
1022	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
1023
1024	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1025	DO l = 1,llm-1
1026	do ij = ijb,ije
1027	IF(w(ij,l+1).gt.0.) THEN
1028	sigw=w(ij,l+1)/masse(ij,l+1)
1029	wq(ij,l+1)=w(ij,l+1)(q(ij,l+1)+0.5(1.-sigw)*dzq(ij,l+1))
1030	ELSE
1031	sigw=w(ij,l+1)/masse(ij,l)
1032	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
1033	ENDIF
1034	ENDDO
1035	ENDDO
1036	c$OMP END DO NOWAIT
1037
1038	c$OMP MASTER
1039	DO ij=ijb,ije
1040	wq(ij,llm+1)=0.
1041	wq(ij,1)=0.
1042	ENDDO
1043	c$OMP END MASTER
1044	c$OMP BARRIER
1045
1046	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1047	DO l=1,llm
1048	DO ij=ijb,ije
1049	newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
1050	q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
1051	& /newmasse
1052	masse(ij,l)=newmasse
1053	ENDDO
1054	ENDDO
1055	c$OMP END DO NOWAIT
1056
1057
1058	RETURN
1059	END
1060	c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
1061	c
1062	c#include "dimensions.h"
1063	c#include "paramet.h"
1064
1065	c CHARACTER() comment
1066	c real qmin,qmax
1067	c real zq(ip1jmp1,llm)
1068
1069	c INTEGER jadrs(ip1jmp1), jbad, k, i
1070
1071
1072	c DO k = 1, llm
1073	c jbad = 0
1074	c DO i = 1, ip1jmp1
1075	c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
1076	c jbad = jbad + 1
1077	c jadrs(jbad) = i
1078	c ENDIF
1079	c ENDDO
1080	c IF (jbad.GT.0) THEN
1081	c PRINT*, comment
1082	c DO i = 1, jbad
1083	cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
1084	c ENDDO
1085	c ENDIF
1086	c ENDDO
1087
1088	c return
1089	c end
1090
1091
1092	subroutine minmaxq_p(zq,qmin,qmax,comment)
1093
1094	#include "dimensions.h"
1095	#include "paramet.h"
1096
1097	character*20 comment
1098	real qmin,qmax
1099	real zq(ip1jmp1,llm)
1100	real zzq(iip1,jjp1,llm)
1101
1102	integer imin,jmin,lmin,ijlmin
1103	integer imax,jmax,lmax,ijlmax
1104
1105	integer ismin,ismax
1106
1107	#ifdef isminismax
1108	call scopy (ip1jmp1*llm,zq,1,zzq,1)
1109
1110	ijlmin=ismin(ijp1llm,zq,1)
1111	lmin=(ijlmin-1)/ip1jmp1+1
1112	ijlmin=ijlmin-(lmin-1.)*ip1jmp1
1113	jmin=(ijlmin-1)/iip1+1
1114	imin=ijlmin-(jmin-1.)*iip1
1115	zqmin=zq(ijlmin,lmin)
1116
1117	ijlmax=ismax(ijp1llm,zq,1)
1118	lmax=(ijlmax-1)/ip1jmp1+1
1119	ijlmax=ijlmax-(lmax-1.)*ip1jmp1
1120	jmax=(ijlmax-1)/iip1+1
1121	imax=ijlmax-(jmax-1.)*iip1
1122	zqmax=zq(ijlmax,lmax)
1123
1124	if(zqmin.lt.qmin)
1125	c s write(*,9999) comment,
1126	s write(,) comment,
1127	s imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
1128	if(zqmax.gt.qmax)
1129	c s write(*,9999) comment,
1130	s write(,) comment,
1131	s imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
1132	#endif
1133	return
1134	9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
1135	end
1136
1137
1138

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

Download in other formats:

Original Format