Context Navigation

vlsplt_loc.F @ 2669

Last change on this file since 2669 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 30.3 KB

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

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

Download in other formats:

Original Format