Context Navigation

vlsplt_loc.F @ 3473

Last change on this file since 3473 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF

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
	1058	c$OMP END DO NO WAIT
	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
	1092	c$OMP END DO NO WAIT
	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