Context Navigation

vlsplt_loc.F @ 4057

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

Fixes for previous commit:

few tracers(:)%nqChilds were errouneously replaced with tarcers(:)%nqDesc
minor changes

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

Rev	Line
[1673]	1	!
	2	! $Id: vlsplt_loc.F 4052 2021-12-26 21:27:36Z dcugnet $
	3	!
[2277]	4	RECURSIVE SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x,iq)
[1632]	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 --------------------------------------------------------------------
[1823]	15	USE parallel_lmdz
[4050]	16	USE infotrac, ONLY : nqtot,tracers, ! CRisi &
[3800]	17	& qperemin,masseqmin,ratiomin ! MVals et CRisi
[1632]	18	IMPLICIT NONE
	19	c
[2597]	20	include "dimensions.h"
	21	include "paramet.h"
[3435]	22	include "iniprint.h"
[1632]	23	c
	24	c
	25	c Arguments:
	26	c ----------
[2270]	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
[1632]	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
[2281]	45	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
[2270]	46	INTEGER ifils,iq2 ! CRisi
	47
[1632]	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
[2286]	57	!write(,) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=',
	58	! & iq,ijb_x
[1632]	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 -----------------------------------------------------
[2270]	72	! on a besoin de q entre ijb et ije
[1632]	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
[2270]	78	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
[1632]	79	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
[2270]	80	c sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
[1632]	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
[2270]	134	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
[1632]	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
[2286]	155	!write(,) 'vlx 156: iq,ijb_x=',iq,ijb_x
[2270]	156
[1632]	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
[2270]	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))
[1632]	181	zdum(ij,l)=0.5*zdum(ij,l)
	182	u_mq(ij,l)=cvmgp(
[2270]	183	, q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
	184	, q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
[1632]	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)
[2270]	195	! on a besoin de masse entre ijb et ije
[1632]	196	DO l=1,llm
	197	DO ij=ijb,ije-1
[2270]	198	c print*,'masse(',ij,')=',masse(ij,l,iq)
[1632]	199	IF (u_m(ij,l).gt.0.) THEN
[2270]	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))
[1632]	203	ELSE
[2270]	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))
[1632]	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'
[2270]	228
[1632]	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'
[2270]	237
[1632]	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)
[2270]	261
	262
[1632]	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
[2286]	274	!PRINT*,'vlx 278, niju,nl',niju,nl(l)
[1632]	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
[2270]	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)
[1632]	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*
[2270]	295	& (q(ijq,l,iq)+0.5*
	296	& (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
[1632]	297	ELSE
	298	ijq=ij+1
	299	i=ijq-(j-1)*iip1
	300	c accumulation pour les mailles completements advectees
[2270]	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)
[1632]	305	i=mod(i,iim)+1
	306	ijq=(j-1)*iip1+i
	307	ENDDO
	308	c ajout de la maille non completement advectee
[2270]	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))
[1632]	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
[2270]	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
[4050]	332	do ifils=1,tracers(iq)%nqDescen
[3800]	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!
[4050]	336	iq2=tracers(iq)%iqDescen(ifils)
[2270]	337	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[4050]	338	DO l=1,llm
[2270]	339	DO ij=ijb,ije
[4050]	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),masseqmin)
	344	if (q(ij,l,iq).gt.qperemin) 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)=ratiomin
	348	endif
[2270]	349	enddo
[4050]	350	enddo
[2270]	351	c$OMP END DO NOWAIT
[4050]	352	enddo !do ifils=1,tracers(iq)%nqDescen
	353	do ifils=1,tracers(iq)%nqChilds
	354	iq2=tracers(iq)%iqDescen(ifils)
	355	call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
	356	enddo
[2270]	357	! end CRisi
	358
	359
[1632]	360	c calcul des tENDances
	361	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	362	DO l=1,llm
	363	DO ij=ijb+1,ije
[3800]	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),masseqmin)
[2270]	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
[1632]	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
[2270]	373	q(ij-iim,l,iq)=q(ij,l,iq)
	374	masse(ij-iim,l,iq)=masse(ij,l,iq)
[1632]	375	ENDDO
	376	ENDDO
	377	c$OMP END DO NOWAIT
[2270]	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
[4050]	382	do ifils=1,tracers(iq)%nqDescen
	383	iq2=tracers(iq)%iqDescen(ifils)
[2270]	384	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[4050]	385	DO l=1,llm
[2270]	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
[4050]	390	q(ij-iim,l,iq2)=q(ij,l,iq2)
	391	enddo
	392	enddo
[2270]	393	c$OMP END DO NOWAIT
[4050]	394	enddo
[2270]	395
[2286]	396	!write(,) 'vlsplt 399: iq,ijb_x=',iq,ijb_x
[1632]	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
[2277]	405	RECURSIVE SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq)
[1632]	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 --------------------------------------------------------------------
[1823]	417	USE parallel_lmdz
[4050]	418	USE infotrac, ONLY : nqtot,tracers, ! CRisi &
[3800]	419	& qperemin,masseqmin,ratiomin ! MVals et CRisi
[2597]	420	USE comconst_mod, ONLY: pi
[1632]	421	IMPLICIT NONE
	422	c
[2597]	423	include "dimensions.h"
	424	include "paramet.h"
	425	include "comgeom.h"
[1632]	426	c
	427	c
	428	c Arguments:
	429	c ----------
[2270]	430	REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
[1632]	431	REAL masse_adv_v( ijb_v:ije_v,llm)
[2270]	432	REAL q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm)
	433	INTEGER iq ! CRisi
[1632]	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
[1673]	445	REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
[1632]	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)
[2270]	462
[2281]	463	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
[2270]	464	INTEGER ifils,iq2 ! CRisi
[1632]	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
[3800]	473	INTEGER ijbm,ijem
[1632]	474
[2270]	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
[1632]	480	IF(first) THEN
[2270]	481	PRINT*,'Shema Amont nouveau appele dans Vanleer '
[1632]	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
[2270]	513	airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
[1632]	514	ENDDO
	515	qpns = SSUM( iim, airescb ,1 ) / airej2
	516	endif
	517
	518	if (pole_sud) then
	519	DO i = 1, iim
[2270]	520	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
[1632]	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
[2270]	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
[1632]	535	DO ij=ijb,ije
[2270]	536	dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
[1632]	537	adyqv(ij)=abs(dyqv(ij))
	538	ENDDO
[2270]	539
[1632]	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
[2270]	556	dyq(ij,l)=qpns-q(ij+iip1,l,iq)
[1632]	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
[2270]	578	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
[1632]	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)
[1673]	627	C appn=abs(dyq(1)/dyqv(iip1+1))
[1632]	628	C PRINT*,dyq(ip1jm+1)
	629	C PRINT*,dyqv(ip1jm-iip1+1)
[1673]	630	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[1632]	631	C DO ij=2,iim
[1673]	632	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
	633	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
[1632]	634	C ENDDO
[1673]	635	C appn=min(pente_max/appn,1.)
	636	C apps=min(pente_max/apps,1.)
[1632]	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.)
[1673]	642	C & appn=0.
[1632]	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.)
[1673]	645	C & apps=0.
[1632]	646	C
	647	C limitation des pentes aux poles
	648	C DO ij=1,iip1
[1673]	649	C dyq(ij)=appn*dyq(ij)
	650	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
[1632]	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
[2270]	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))
[1632]	717	ELSE
[2270]	718	qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
	719	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq))
[1632]	720	ENDIF
	721	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
	722	ENDDO
	723	ENDDO
	724	c$OMP END DO NOWAIT
[2270]	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
[4052]	728	!write(,) 'vly 689: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
[2270]	729
	730	ijb=ij_begin-2*iip1
	731	ije=ij_end+2*iip1
[3800]	732	ijbm=ij_begin-iip1
	733	ijem=ij_end+iip1
[2270]	734	if (pole_nord) ijb=ij_begin
[3800]	735	if (pole_sud) ije=ij_end
	736	if (pole_nord) ijbm=ij_begin
	737	if (pole_sud) ijem=ij_end
	738
[4050]	739	do ifils=1,tracers(iq)%nqDescen
	740	iq2=tracers(iq)%iqDescen(ifils)
[2270]	741	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[4050]	742	DO l=1,llm
	743	! modif des bornes: CRisi 16 nov 2020
	744	! d'abord masse avec bornes corrigées
[3800]	745	DO ij=ijbm,ijem
[4050]	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),masseqmin)
	748	enddo
[3800]	749
	750	! ensuite Ratio avec anciennes bornes
[4050]	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.qperemin) 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)=ratiomin
	757	endif
[3800]	758	enddo !DO ij=ijbm,ijem
[4050]	759	enddo !DO l=1,llm
[2270]	760	c$OMP END DO NOWAIT
[4050]	761	enddo
[2270]	762
[4050]	763	do ifils=1,tracers(iq)%nqChilds
	764	iq2=tracers(iq)%iqDescen(ifils)
	765	call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
	766	enddo
[2270]	767	! end CRisi
[1632]	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
[2270]	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
[1632]	785	ENDDO
[2270]	786
	787
[1632]	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)
[2270]	794	massepn=ssum(iim,masse(1,l,iq),1)
[1632]	795	qpn=0.
	796	do ij=1,iim
[2270]	797	qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
[1632]	798	enddo
	799	qpn=(qpn+convpn)/(massepn+convmpn)
	800	do ij=1,iip1
[2270]	801	q(ij,l,iq)=qpn
[1632]	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)
[2270]	812	masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
[1632]	813	qps=0.
	814	do ij = ip1jm+1,ip1jmp1-1
[2270]	815	qps=qps+masse(ij,l,iq)*q(ij,l,iq)
[1632]	816	enddo
	817	qps=(qps+convps)/(masseps+convmps)
	818	do ij=ip1jm+1,ip1jmp1
[2270]	819	q(ij,l,iq)=qps
[1632]	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
[2270]	833	c masse(ij,l,iq)=newmasse*aire(ij)
[1632]	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
[2270]	843	c masse(ij,l,iq)=newmasse*aire(ij)
[1632]	844	c ENDDO
	845	c._. fin nouvelle version
	846	ENDDO
	847	c$OMP END DO NOWAIT
	848
[2270]	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
[4050]	855	do ifils=1,tracers(iq)%nqDescen
	856	iq2=tracers(iq)%iqDescen(ifils)
[2270]	857	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[4050]	858	DO l=1,llm
[2270]	859	DO ij=ijb,ije
	860	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
	861	enddo
[4050]	862	enddo
[2270]	863	c$OMP END DO NOWAIT
[4050]	864	enddo
[2270]	865
	866
[1632]	867	RETURN
	868	END
	869
	870
	871
[2277]	872	RECURSIVE SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq)
[1632]	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 --------------------------------------------------------------------
[1823]	884	USE parallel_lmdz
[1632]	885	USE vlz_mod
[4050]	886	USE infotrac, ONLY : nqtot,tracers, ! CRisi &
[3800]	887	& qperemin,masseqmin,ratiomin ! MVals et CRisi
	888
[1632]	889	IMPLICIT NONE
	890	c
[2597]	891	include "dimensions.h"
	892	include "paramet.h"
[3435]	893	include "iniprint.h"
[1632]	894	c
	895	c
	896	c Arguments:
	897	c ----------
[2270]	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
[1632]	902	c
	903	c Local
	904	c ---------
	905	c
	906	INTEGER i,ij,l,j,ii
[2765]	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)
[1632]	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
[2270]	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
[2765]	939
[1632]	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
[2286]	946	!write(,) 'vlsplt 926: entree dans vlz_loc, iq=',iq
[1632]	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
[2270]	959	dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
[1632]	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
[2765]	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
[1632]	1013	c ---------------------------------------------------------------
[2765]	1014	c Identification des mailles ou on viole le CFL : w > masse
	1015	c ---------------------------------------------------------------
	1016
	1017	IF (countcfl==0) THEN
	1018
	1019	c ---------------------------------------------------------------
[1632]	1020	c .... calcul des termes d'advection verticale .......
[2765]	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
[1632]	1025	c ---------------------------------------------------------------
	1026
	1027	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
	1028
[2286]	1029	!write(,) 'vlz 982,ijb,ije=',ijb,ije
[1632]	1030	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	1031	DO l = 1,llm-1
	1032	do ij = ijb,ije
[2270]	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))
[1632]	1037	ELSE
[2270]	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))
[1632]	1041	ENDIF
	1042	ENDDO
	1043	ENDDO
[2270]	1044	c$OMP END DO NOWAIT
[2286]	1045	!write(,) 'vlz 1001'
[1632]	1046
[2765]	1047	ELSE ! countcfl>=1
	1048
[3435]	1049	IF (prt_level>9) THEN
	1050	WRITE(lunout,*)'vlz passage dans le non local'
	1051	ENDIF
[2765]	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 NO WAIT
	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)
[3435]	1083	IF (prt_level>9) THEN
	1084	WRITE(lunout,*)'On viole le CFL Vertical sur ',countcfl,' pts'
	1085	ENDIF
[2765]	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
[3800]	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
[2765]	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 NO WAIT
	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
[1632]	1140	c$OMP MASTER
	1141	DO ij=ijb,ije
[2270]	1142	wq(ij,llm+1,iq)=0.
	1143	wq(ij,1,iq)=0.
[1632]	1144	ENDDO
	1145	c$OMP END MASTER
	1146	c$OMP BARRIER
	1147
[2270]	1148	! CRisi: appel récursif de l'advection sur les fils.
	1149	! Il faut faire ça avant d'avoir mis à jour q et masse
[4050]	1150	!write(,) 'vlsplt 942: iq,nqChilds(iq)=',iq,tracers(iq)%nqChilds
	1151	do ifils=1,tracers(iq)%nqDescen
	1152	iq2=tracers(iq)%iqDescen(ifils)
[1632]	1153	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[4050]	1154	DO l=1,llm
[2270]	1155	DO ij=ijb,ije
[3800]	1156	!MVals: veiller a ce qu'on n'ait pas de denominateur nul
[4050]	1157	masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),masseqmin)
	1158	if (q(ij,l,iq).gt.qperemin) then
	1159	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
	1160	else
	1161	Ratio(ij,l,iq2)=ratiomin
	1162	endif
	1163	!wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
	1164	w(ij,l,iq2)=wq(ij,l,iq)
[2270]	1165	enddo
[4050]	1166	enddo
[2270]	1167	c$OMP END DO NOWAIT
[4050]	1168	enddo
[2270]	1169	c$OMP BARRIER
[2281]	1170
[4050]	1171	do ifils=1,tracers(iq)%nqChilds
	1172	iq2=tracers(iq)%iqDescen(ifils)
	1173	call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2)
	1174	enddo
[2270]	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
[2281]	1179
[2270]	1180	c$OMP BARRIER
	1181	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1182	DO l=1,llm
	1183	DO ij=ijb,ije
[2270]	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))
[1632]	1187	& /newmasse
[2270]	1188	masse(ij,l,iq)=newmasse
[1632]	1189	ENDDO
	1190	ENDDO
	1191	c$OMP END DO NOWAIT
	1192
[2281]	1193
[2270]	1194	! retablir les fils en rapport de melange par rapport a l'air:
[4050]	1195	do ifils=1,tracers(iq)%nqDescen
	1196	iq2=tracers(iq)%iqDescen(ifils)
[2270]	1197	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[4050]	1198	DO l=1,llm
[2270]	1199	DO ij=ijb,ije
	1200	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
	1201	enddo
[4050]	1202	enddo
[2270]	1203	c$OMP END DO NOWAIT
[4050]	1204	enddo
[1632]	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