Context Navigation

vlsplt_p.F @ 1824

Last change on this file since 1824 was 1823, checked in by Ehouarn Millour, 11 years ago
Remplacement de parallel.F90 (en conflit avec orchidée) par parallel_lmdz.F90. UG ......................................... Renaming parallel.F90 (conflicting with orchidée) into parallel_lmdz.F90. UG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.9 KB

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

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

Download in other formats:

Original Format