Context Navigation

vlsplt_p.F @ 5223

Last change on this file since 5223 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 29.9 KB

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

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

Download in other formats:

Original Format