Context Navigation

vlsplt_loc.F @ 4607

Last change on this file since 4607 was 4593, checked in by yann meurdesoif, 17 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM

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

File size: 35.2 KB

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

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

Download in other formats:

Original Format