Context Navigation

vlsplt_loc.F @ 3852

Last change on this file since 3852 was 3852, checked in by dcugnet, 3 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

the number of tracers on the first line
one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

Format 2: read the sections from the "tracer.def" file.
format 3: read one section each "tracer_<section name>.def" file.
the first line of a section is "&<section name>
the other lines start with a tracer name followed by <key>=<val> pairs.
the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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.9 KB

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

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

Download in other formats:

Original Format