Context Navigation

vlsplt.F @ 3999

Last change on this file since 3999 was 3852, checked in by dcugnet, 4 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

the number of tracers on the first line
one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

Format 2: read the sections from the "tracer.def" file.
format 3: read one section each "tracer_<section name>.def" file.
the first line of a section is "&<section name>
the other lines start with a tracer name followed by <key>=<val> pairs.
the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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

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

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

Download in other formats:

Original Format