Context Navigation

vlsplt.F @ 3670

Last change on this file since 3670 was 357, checked in by lmdz, 23 years ago
Remplacement, dans certaines circonstances l'ancienne version pouvait faire planter le modele. FH LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 25.1 KB

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

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

Download in other formats:

Original Format