Context Navigation

vlsplt.F @ 1573

Last change on this file since 1573 was 1550, checked in by lguez, 13 years ago

Bug fix in "bilan_dyn_p". The index was out of bounds in the removed
assignment . Also, the removed assignment was useless.

Bug fix in "coefkzmin". The size of a dummy array cannot exceed the
size of the associated actual array. ("coefkzmin" is called by
"coef_diff_turb".) "km(:, klev+1)" and "kn(:, klev+1)" were not
defined in "coefkzmin" so this was maybe an innocuous bug.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 25.2 KB

Rev	Line
[524]	1	c
[1520]	2	c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z jghattas $
[524]	3	c
	4
	5	SUBROUTINE vlsplt(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 --------------------------------------------------------------------
	20	IMPLICIT NONE
	21	c
	22	#include "dimensions.h"
	23	#include "paramet.h"
	24	#include "logic.h"
	25	#include "comvert.h"
	26	#include "comconst.h"
	27
	28	c
	29	c Arguments:
	30	c ----------
	31	REAL masse(ip1jmp1,llm),pente_max
	32	c REAL masse(iip1,jjp1,llm),pente_max
	33	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
	34	REAL q(ip1jmp1,llm)
	35	c REAL q(iip1,jjp1,llm)
	36	REAL w(ip1jmp1,llm),pdt
	37	c
	38	c Local
	39	c ---------
	40	c
	41	INTEGER i,ij,l,j,ii
	42	INTEGER ijlqmin,iqmin,jqmin,lqmin
	43	c
	44	REAL zm(ip1jmp1,llm),newmasse
	45	REAL mu(ip1jmp1,llm)
	46	REAL mv(ip1jm,llm)
	47	REAL mw(ip1jmp1,llm+1)
	48	REAL zq(ip1jmp1,llm),zz
	49	REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
	50	REAL second,temps0,temps1,temps2,temps3
	51	REAL ztemps1,ztemps2,ztemps3
	52	REAL zzpbar, zzw
	53	LOGICAL testcpu
	54	SAVE testcpu
	55	SAVE temps1,temps2,temps3
	56	INTEGER iminn,imaxx
	57
	58	REAL qmin,qmax
	59	DATA qmin,qmax/0.,1.e33/
	60	DATA testcpu/.false./
	61	DATA temps1,temps2,temps3/0.,0.,0./
	62
	63
	64	zzpbar = 0.5 * pdt
	65	zzw = pdt
	66	DO l=1,llm
	67	DO ij = iip2,ip1jm
	68	mu(ij,l)=pbaru(ij,l) * zzpbar
	69	ENDDO
	70	DO ij=1,ip1jm
	71	mv(ij,l)=pbarv(ij,l) * zzpbar
	72	ENDDO
	73	DO ij=1,ip1jmp1
	74	mw(ij,l)=w(ij,l) * zzw
	75	ENDDO
	76	ENDDO
	77
	78	DO ij=1,ip1jmp1
	79	mw(ij,llm+1)=0.
	80	ENDDO
	81
	82	CALL SCOPY(ijp1llm,q,1,zq,1)
	83	CALL SCOPY(ijp1llm,masse,1,zm,1)
	84
	85	cprint*,'Entree vlx1'
	86	c call minmaxq(zq,qmin,qmax,'avant vlx ')
	87	call vlx(zq,pente_max,zm,mu)
	88	cprint*,'Sortie vlx1'
	89	c call minmaxq(zq,qmin,qmax,'apres vlx1 ')
	90
	91	c print*,'Entree vly1'
	92	call vly(zq,pente_max,zm,mv)
	93	c call minmaxq(zq,qmin,qmax,'apres vly1 ')
	94	cprint*,'Sortie vly1'
	95	call vlz(zq,pente_max,zm,mw)
	96	c call minmaxq(zq,qmin,qmax,'apres vlz ')
	97
	98
	99	call vly(zq,pente_max,zm,mv)
	100	c call minmaxq(zq,qmin,qmax,'apres vly ')
	101
	102
	103	call vlx(zq,pente_max,zm,mu)
	104	c call minmaxq(zq,qmin,qmax,'apres vlx2 ')
	105
	106
	107	DO l=1,llm
	108	DO ij=1,ip1jmp1
	109	q(ij,l)=zq(ij,l)
	110	ENDDO
	111	DO ij=1,ip1jm+1,iip1
	112	q(ij+iim,l)=q(ij,l)
	113	ENDDO
	114	ENDDO
	115
	116	RETURN
	117	END
	118	SUBROUTINE vlx(q,pente_max,masse,u_m)
	119
	120	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	121	c
	122	c ********************************************************************
	123	c Shema d'advection " pseudo amont " .
	124	c ********************************************************************
	125	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	126	c
	127	c
	128	c --------------------------------------------------------------------
	129	IMPLICIT NONE
	130	c
[1550]	131	include "dimensions.h"
	132	include "paramet.h"
	133	include "logic.h"
	134	include "comvert.h"
	135	include "comconst.h"
	136	include "iniprint.h"
[524]	137	c
	138	c
	139	c Arguments:
	140	c ----------
	141	REAL masse(ip1jmp1,llm),pente_max
	142	REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
	143	REAL q(ip1jmp1,llm)
	144	REAL w(ip1jmp1,llm)
	145	c
	146	c Local
	147	c ---------
	148	c
	149	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
	150	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
	151	c
	152	REAL new_m,zu_m,zdum(ip1jmp1,llm)
	153	REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
	154	REAL zz(ip1jmp1)
	155	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
	156	REAL u_mq(ip1jmp1,llm)
	157
	158	Logical extremum,first,testcpu
	159	SAVE first,testcpu
	160
	161	REAL SSUM
	162	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
	163	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	164
	165	REAL z1,z2,z3
	166
	167	DATA first,testcpu/.true.,.false./
	168
	169	IF(first) THEN
	170	temps1=0.
	171	temps2=0.
	172	temps3=0.
	173	temps4=0.
	174	temps5=0.
	175	first=.false.
	176	ENDIF
	177
	178	c calcul de la pente a droite et a gauche de la maille
	179
	180
	181	IF (pente_max.gt.-1.e-5) THEN
	182	c IF (pente_max.gt.10) THEN
	183
	184	c calcul des pentes avec limitation, Van Leer scheme I:
	185	c -----------------------------------------------------
	186
	187	c calcul de la pente aux points u
	188	DO l = 1, llm
	189	DO ij=iip2,ip1jm-1
	190	dxqu(ij)=q(ij+1,l)-q(ij,l)
	191	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	192	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	193	ENDDO
	194	DO ij=iip1+iip1,ip1jm,iip1
	195	dxqu(ij)=dxqu(ij-iim)
	196	c sigu(ij)=sigu(ij-iim)
	197	ENDDO
	198
	199	DO ij=iip2,ip1jm
	200	adxqu(ij)=abs(dxqu(ij))
	201	ENDDO
	202
	203	c calcul de la pente maximum dans la maille en valeur absolue
	204
	205	DO ij=iip2+1,ip1jm
	206	dxqmax(ij,l)=pente_max*
	207	, min(adxqu(ij-1),adxqu(ij))
	208	c limitation subtile
	209	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	210
	211
	212	ENDDO
	213
	214	DO ij=iip1+iip1,ip1jm,iip1
	215	dxqmax(ij-iim,l)=dxqmax(ij,l)
	216	ENDDO
	217
	218	DO ij=iip2+1,ip1jm
	219	#ifdef CRAY
	220	dxq(ij,l)=
	221	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	222	#else
	223	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	224	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	225	ELSE
	226	c extremum local
	227	dxq(ij,l)=0.
	228	ENDIF
	229	#endif
	230	dxq(ij,l)=0.5*dxq(ij,l)
	231	dxq(ij,l)=
	232	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	233	ENDDO
	234
	235	ENDDO ! l=1,llm
	236	cprint*,'Ok calcul des pentes'
	237
	238	ELSE ! (pente_max.lt.-1.e-5)
	239
	240	c Pentes produits:
	241	c ----------------
	242
	243	DO l = 1, llm
	244	DO ij=iip2,ip1jm-1
	245	dxqu(ij)=q(ij+1,l)-q(ij,l)
	246	ENDDO
	247	DO ij=iip1+iip1,ip1jm,iip1
	248	dxqu(ij)=dxqu(ij-iim)
	249	ENDDO
	250
	251	DO ij=iip2+1,ip1jm
	252	zz(ij)=dxqu(ij-1)*dxqu(ij)
	253	zz(ij)=zz(ij)+zz(ij)
	254	IF(zz(ij).gt.0) THEN
	255	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	256	ELSE
	257	c extremum local
	258	dxq(ij,l)=0.
	259	ENDIF
	260	ENDDO
	261
	262	ENDDO
	263
	264	ENDIF ! (pente_max.lt.-1.e-5)
	265
	266	c bouclage de la pente en iip1:
	267	c -----------------------------
	268
	269	DO l=1,llm
	270	DO ij=iip1+iip1,ip1jm,iip1
	271	dxq(ij-iim,l)=dxq(ij,l)
	272	ENDDO
	273	DO ij=1,ip1jmp1
	274	iadvplus(ij,l)=0
	275	ENDDO
	276
	277	ENDDO
	278
	279	c print*,'Bouclage en iip1'
	280
	281	c calcul des flux a gauche et a droite
	282
	283	#ifdef CRAY
	284
	285	DO l=1,llm
	286	DO ij=iip2,ip1jm-1
	287	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	288	, 1.+u_m(ij,l)/masse(ij+1,l),
	289	, u_m(ij,l))
	290	zdum(ij,l)=0.5*zdum(ij,l)
	291	u_mq(ij,l)=cvmgp(
	292	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	293	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	294	, u_m(ij,l))
	295	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	296	ENDDO
	297	ENDDO
	298	#else
	299	c on cumule le flux correspondant a toutes les mailles dont la masse
	300	c au travers de la paroi pENDant le pas de temps.
	301	cprint*,'Cumule ....'
	302
	303	DO l=1,llm
	304	DO ij=iip2,ip1jm-1
	305	c print*,'masse(',ij,')=',masse(ij,l)
	306	IF (u_m(ij,l).gt.0.) THEN
	307	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	308	u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
	309	ELSE
	310	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	311	u_mq(ij,l)=u_m(ij,l)(q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
	312	ENDIF
	313	ENDDO
	314	ENDDO
	315	#endif
	316	c stop
	317
	318	c go to 9999
	319	c detection des points ou on advecte plus que la masse de la
	320	c maille
	321	DO l=1,llm
	322	DO ij=iip2,ip1jm-1
	323	IF(zdum(ij,l).lt.0) THEN
	324	iadvplus(ij,l)=1
	325	u_mq(ij,l)=0.
	326	ENDIF
	327	ENDDO
	328	ENDDO
	329	cprint*,'Ok test 1'
	330	DO l=1,llm
	331	DO ij=iip1+iip1,ip1jm,iip1
	332	iadvplus(ij,l)=iadvplus(ij-iim,l)
	333	ENDDO
	334	ENDDO
	335	c print*,'Ok test 2'
	336
	337
	338	c traitement special pour le cas ou on advecte en longitude plus que le
	339	c contenu de la maille.
	340	c cette partie est mal vectorisee.
	341
	342	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	343
	344	n0=0
	345	DO l=1,llm
	346	nl(l)=0
	347	DO ij=iip2,ip1jm
	348	nl(l)=nl(l)+iadvplus(ij,l)
	349	ENDDO
	350	n0=n0+nl(l)
	351	ENDDO
	352
[595]	353	IF(n0.gt.0) THEN
[1550]	354	if (prt_level > 2) PRINT *,
	355	$ 'Nombre de points pour lesquels on advect plus que le'
[524]	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
[1520]	480	REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
[524]	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)
[1520]	604	C appn=abs(dyq(1)/dyqv(iip1+1))
[524]	605	C PRINT*,dyq(ip1jm+1)
	606	C PRINT*,dyqv(ip1jm-iip1+1)
[1520]	607	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[524]	608	C DO ij=2,iim
[1520]	609	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
	610	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
[524]	611	C ENDDO
[1520]	612	C appn=min(pente_max/appn,1.)
	613	C apps=min(pente_max/apps,1.)
[524]	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.)
[1520]	619	C & appn=0.
[524]	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.)
[1520]	622	C & apps=0.
[524]	623	C
	624	C limitation des pentes aux poles
	625	C DO ij=1,iip1
[1520]	626	C dyq(ij)=appn*dyq(ij)
	627	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
[524]	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