Context Navigation

vlsplt.F @ 4113

Last change on this file since 4113 was 4064, checked in by dcugnet, 3 years ago

minor fixes (unused variables suppressed, comas after a WRITE() statement, etc.)
parser routines taken from version 7 of https://svn.lmd.jussieu.fr/tracers-parser
few changes in infotrac, and few fixes of (at least) the sequential version:
- uadv and vadv were deallocated twice (fix was lost by mistake just before last commit)
- in [( dum(im), im=1, nm)] implicit loops, ifort evaluates "dum(im)" even if nm==0,

resulting in a crash, "im" being unitialized.

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.1 KB

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

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

Download in other formats:

Original Format