Context Navigation

vlsplt.F @ 2919

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

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

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

Download in other formats:

Original Format