Context Navigation

vlsplt_p.F @ 3000

Last change on this file since 3000 was 2307, checked in by mvals, 5 years ago
Mars GCM: follow-up of the commit regarding the dynamical transport of isotopes: new variables for the thresholds for zq(pere) and masseq in the transport of the isotopic Ratio MV
File size: 37.5 KB

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

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

Download in other formats:

Original Format