Context Navigation

vlsplt_loc.F @ 5049

Last change on this file since 5049 was 3703, checked in by adurocher, 4 years ago

Fixed compilation errors

Fixed some OpenMP syntax errors for scorep
Fix issues with ifort -check all -warn all -O0

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

File size: 34.0 KB

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

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

Download in other formats:

Original Format