Context Navigation

vlsplt.F @ 3492

Last change on this file since 3492 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 30.6 KB

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

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

Download in other formats:

Original Format