Context Navigation

vlsplt_loc.F @ 2410

Last change on this file since 2410 was 2286, checked in by Ehouarn Millour, 10 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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

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

Download in other formats:

Original Format