Context Navigation

vlsplt_loc.F @ 4412

Last change on this file since 4412 was 4325, checked in by dcugnet, 2 years ago
simplify the parser usage: the getKey_init routine is now embedded in the readTracersFile routine. the initIsotopes routine is now embedded in the readIsotopesFile routine. the database is now unique, but can be changed using the get/setKeysDBase. the derived types descriptions, originally located in trac_types_mod, are moved to readTracFiles_mod. few checkings moved from infotrac to the routine testIsotopes, contained in the readIsotopesFile function from readTracFiles_mod. the readTracersFiles and readIsotopesFile routines no longer use a tracers/isotopes argument. remove tnat and alpha_ideal from infotrac ; use instead getKey to get them where they are used (check_isotopes, dynetat0, iniacademic) the trac_type field %Childs is renamed %Children move the isoSelect routine and the corresponding variables routine from infotrac and infotrac_phy to readTracFiles_mod infotrac_phy routine is now fully independant of the (very similar) routine infotrac (init_infotrac_phy has no arguments left). all the explicit keys of the trac_type are now included in the embedded keys database, accessible using the getKey function. the getKey/addKey routines are expanded to handle vectors of integers, reals, logicals or strings. few subroutines converted into functions with error return value. corrections for isotopic tagging tracers mode (to be continued).
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:keywords set to `Id`
File size: 35.2 KB

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

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

Download in other formats:

Original Format