Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

vlsplt_p.F

Last change on this file was 2603, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn logic.h into module logic_mod.F90 EM
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.6 KB

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

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

Download in other formats:

Original Format