Context Navigation

vlsplt_loc.F @ 3687

Last change on this file since 3687 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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

File size: 34.0 KB

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

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

Download in other formats:

Original Format