Context Navigation

vlsplt_p.F @ 1443

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

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

Download in other formats:

Original Format