Context Navigation

vlsplt_loc.F @ 2396

Last change on this file since 2396 was 2298, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes -r2237:2291 into testing branch
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
[1707]	1	!
	2	! $Id$
	3	!
[2298]	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 --------------------------------------------------------------------
[1864]	15	USE parallel_lmdz
[2298]	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 ----------
[2298]	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
[2298]	46	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
	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
[2298]	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 -----------------------------------------------------
[2298]	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
[2298]	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'
[2298]	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
[2298]	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
[2298]	156	!write(,) 'vlx 156: iq,ijb_x=',iq,ijb_x
	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
[2298]	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(
[2298]	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)
[2298]	196	! on a besoin de masse entre ijb et ije
[1632]	197	DO l=1,llm
	198	DO ij=ijb,ije-1
[2298]	199	c print*,'masse(',ij,')=',masse(ij,l,iq)
[1632]	200	IF (u_m(ij,l).gt.0.) THEN
[2298]	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
[2298]	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'
[2298]	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'
[2298]	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)
[2298]	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
[2298]	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
[2298]	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*
[2298]	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
[2298]	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
[2298]	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
[2298]	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
	333	!write(,) 'vlsplt 326: iq,ijb_x,nqfils(iq)=',iq,ijb_x,nqfils(iq)
	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
	343	masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
	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)
	351	call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
	352	enddo !do ifils=1,nqfils(iq)
	353	endif !if (nqfils(iq).gt.0) then
	354	! end CRisi
	355
	356	!write(,) 'vlsplt 360: iq,ijb_x=',iq,ijb_x
	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
[2298]	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
[2298]	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
[2298]	375	!write(,) 'vlsplt 380: iq,ijb_x=',iq,ijb_x
	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
	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
[2298]	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 --------------------------------------------------------------------
[1864]	417	USE parallel_lmdz
[2298]	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 ----------
[2298]	431	REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
[1632]	432	REAL masse_adv_v( ijb_v:ije_v,llm)
[2298]	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
[1707]	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)
[2298]	463
	464	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	536	dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
[1632]	537	adyqv(ij)=abs(dyqv(ij))
	538	ENDDO
[2298]	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
[2298]	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
[2298]	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)
[1707]	627	C appn=abs(dyq(1)/dyqv(iip1+1))
[1632]	628	C PRINT*,dyq(ip1jm+1)
	629	C PRINT*,dyqv(ip1jm-iip1+1)
[1707]	630	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[1632]	631	C DO ij=2,iim
[1707]	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
[1707]	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.)
[1707]	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.)
[1707]	645	C & apps=0.
[1632]	646	C
	647	C limitation des pentes aux poles
	648	C DO ij=1,iip1
[1707]	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
[2298]	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
[2298]	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
[2298]	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
	728	!write(,) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
	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
	741	masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
	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)
	750	call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
	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
[2298]	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
[2298]	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)
[2298]	780	massepn=ssum(iim,masse(1,l,iq),1)
[1632]	781	qpn=0.
	782	do ij=1,iim
[2298]	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
[2298]	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)
[2298]	798	masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
[1632]	799	qps=0.
	800	do ij = ip1jm+1,ip1jmp1-1
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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 --------------------------------------------------------------------
[1864]	872	USE parallel_lmdz
[1632]	873	USE vlz_mod
[2298]	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 ----------
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	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
[2298]	999	c$OMP END DO NOWAIT
	1000	!write(,) 'vlz 1001'
[1632]	1001
	1002	c$OMP MASTER
	1003	DO ij=ijb,ije
[2298]	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
[2298]	1010	! CRisi: appel récursif de l'advection sur les fils.
	1011	! Il faut faire ça avant d'avoir mis à jour q et masse
	1012	!write(,) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
	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)
[2298]	1017	DO l=1,llm
	1018	DO ij=ijb,ije
	1019	masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
	1020	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
	1021	!wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
	1022	w(ij,l,iq2)=wq(ij,l,iq)
	1023	enddo
	1024	enddo
	1025	c$OMP END DO NOWAIT
	1026	enddo !do ifils=1,nqdesc(iq)
	1027	c$OMP BARRIER
	1028
	1029	do ifils=1,nqfils(iq)
	1030	iq2=iqfils(ifils,iq)
	1031	call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2)
	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
	1038
	1039	c$OMP BARRIER
	1040	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1632]	1041	DO l=1,llm
	1042	DO ij=ijb,ije
[2298]	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
[2298]	1047	masse(ij,l,iq)=newmasse
[1632]	1048	ENDDO
	1049	ENDDO
	1050	c$OMP END DO NOWAIT
	1051
[2298]	1052
	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