Context Navigation

vlsplt.F @ 3173

Last change on this file since 3173 was 595, checked in by Laurent Fairhead, 20 years ago
Si il n'y avait qu'un point ou on advectait plus que le contenu de la maille, il n'etait jamais traite YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.1 KB

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

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

Download in other formats:

Original Format