Context Navigation

vlsplt.F @ 5006

Last change on this file since 5006 was 4470, checked in by Laurent Fairhead, 22 months ago
Replaced STOP instructions by calls to abort_gcm for a cleaner exit
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.0 KB

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

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

Download in other formats:

Original Format