Context Navigation

vlsplt_loc.F @ 5452

Last change on this file since 5452 was 2765, checked in by fhourdin, 8 years ago

Inclusion d'un traitement spécifique dans l'advection verticale
pour les cas de nombre de courant w/masse > 1
ou w est le transfert de masse entre deux couches.

ce traitement garanti la convergence numérique quand le critère
n'est jamais violé.
Dans les autres cas, il fait proprement le transport vertical
qui est autrement hors de contrôle.
Pourrait contribuer a resoudre les plantages himalayens.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 33.8 KB

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

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

Download in other formats:

Original Format