Context Navigation

vlsplt_p.F @ 2902

Last change on this file since 2902 was 2307, checked in by mvals, 5 years ago
Mars GCM: follow-up of the commit regarding the dynamical transport of isotopes: new variables for the thresholds for zq(pere) and masseq in the transport of the isotopic Ratio MV
File size: 37.5 KB

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

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

Download in other formats:

Original Format