Context Navigation

vlspltqs.F @ 5305

Last change on this file since 5305 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: 21.1 KB

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
	5	, p,pk,teta )
	6	c
	7	c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
	8	c
	9	c ********************************************************************
	10	c Shema d'advection " pseudo amont " .
	11	c + test sur humidite specifique: Q advecte< Qsat aval
	12	c (F. Codron, 10/99)
	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 teta temperature potentielle, p pression aux interfaces,
	22	c pk exner au milieu des couches necessaire pour calculer Qsat
	23	c --------------------------------------------------------------------
	24	IMPLICIT NONE
	25	c
	26	#include "dimensions.h"
	27	#include "paramet.h"
	28	#include "logic.h"
	29	#include "comvert.h"
	30	#include "comconst.h"
	31
	32	c
	33	c Arguments:
	34	c ----------
	35	REAL masse(ip1jmp1,llm),pente_max
	36	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
	37	REAL q(ip1jmp1,llm)
	38	REAL w(ip1jmp1,llm),pdt
	39	REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
	40	c
	41	c Local
	42	c ---------
	43	c
	44	INTEGER i,ij,l,j,ii
	45	c
	46	REAL qsat(ip1jmp1,llm)
	47	REAL zm(ip1jmp1,llm)
	48	REAL mu(ip1jmp1,llm)
	49	REAL mv(ip1jm,llm)
	50	REAL mw(ip1jmp1,llm+1)
	51	REAL zq(ip1jmp1,llm)
	52	REAL temps1,temps2,temps3
	53	REAL zzpbar, zzw
	54	LOGICAL testcpu
	55	SAVE testcpu
	56	SAVE temps1,temps2,temps3
	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	c--pour rapport de melange saturant--
	64
	65	REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
	66	REAL ptarg,pdelarg,foeew,zdelta
	67	REAL tempe(ip1jmp1)
	68
	69	c fonction psat(T)
	70
	71	FOEEW ( PTARG,PDELARG ) = EXP (
	72	* (R3LES(1.-PDELARG)+R3IESPDELARG) * (PTARG-RTT)
	73	* / (PTARG-(R4LES(1.-PDELARG)+R4IESPDELARG)) )
	74
	75	r2es = 380.11733
	76	r3les = 17.269
	77	r3ies = 21.875
	78	r4les = 35.86
	79	r4ies = 7.66
	80	retv = 0.6077667
	81	rtt = 273.16
	82
	83	c-- Calcul de Qsat en chaque point
	84	c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
	85	c pour eviter une exponentielle.
	86	DO l = 1, llm
	87	DO ij = 1, ip1jmp1
	88	tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
	89	ENDDO
	90	DO ij = 1, ip1jmp1
	91	zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
	92	play = 0.5*(p(ij,l)+p(ij,l+1))
	93	qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
	94	qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
	95	ENDDO
	96	ENDDO
	97
	98	c PRINT*,'Debut vlsplt version debug sans vlyqs'
	99
	100	zzpbar = 0.5 * pdt
	101	zzw = pdt
	102	DO l=1,llm
	103	DO ij = iip2,ip1jm
	104	mu(ij,l)=pbaru(ij,l) * zzpbar
	105	ENDDO
	106	DO ij=1,ip1jm
	107	mv(ij,l)=pbarv(ij,l) * zzpbar
	108	ENDDO
	109	DO ij=1,ip1jmp1
	110	mw(ij,l)=w(ij,l) * zzw
	111	ENDDO
	112	ENDDO
	113
	114	DO ij=1,ip1jmp1
	115	mw(ij,llm+1)=0.
	116	ENDDO
	117
	118	CALL SCOPY(ijp1llm,q,1,zq,1)
	119	CALL SCOPY(ijp1llm,masse,1,zm,1)
	120
	121	c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
	122	call vlxqs(zq,pente_max,zm,mu,qsat)
	123
	124
	125	c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
	126
	127	call vlyqs(zq,pente_max,zm,mv,qsat)
	128
	129
	130	c call minmaxq(zq,qmin,qmax,'avant vlz ')
	131
	132	call vlz(zq,pente_max,zm,mw)
	133
	134
	135	c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
	136	c call minmaxq(zm,qmin,qmax,'M avant vlyqs ')
	137
	138	call vlyqs(zq,pente_max,zm,mv,qsat)
	139
	140
	141	c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
	142	c call minmaxq(zm,qmin,qmax,'M avant vlxqs ')
	143
	144	call vlxqs(zq,pente_max,zm,mu,qsat)
	145
	146	c call minmaxq(zq,qmin,qmax,'apres vlxqs ')
	147	c call minmaxq(zm,qmin,qmax,'M apres vlxqs ')
	148
	149
	150	DO l=1,llm
	151	DO ij=1,ip1jmp1
	152	q(ij,l)=zq(ij,l)
	153	ENDDO
	154	DO ij=1,ip1jm+1,iip1
	155	q(ij+iim,l)=q(ij,l)
	156	ENDDO
	157	ENDDO
	158
	159	RETURN
	160	END
	161	SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
	162	c
	163	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	164	c
	165	c ********************************************************************
	166	c Shema d'advection " pseudo amont " .
	167	c ********************************************************************
	168	c
	169	c --------------------------------------------------------------------
	170	IMPLICIT NONE
	171	c
	172	#include "dimensions.h"
	173	#include "paramet.h"
	174	#include "logic.h"
	175	#include "comvert.h"
	176	#include "comconst.h"
	177	c
	178	c
	179	c Arguments:
	180	c ----------
	181	REAL masse(ip1jmp1,llm),pente_max
	182	REAL u_m( ip1jmp1,llm )
	183	REAL q(ip1jmp1,llm)
	184	REAL qsat(ip1jmp1,llm)
	185	c
	186	c Local
	187	c ---------
	188	c
	189	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
	190	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
	191	c
	192	REAL new_m,zu_m,zdum(ip1jmp1,llm)
	193	REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
	194	REAL zz(ip1jmp1)
	195	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
	196	REAL u_mq(ip1jmp1,llm)
	197
	198	Logical first,testcpu
	199	SAVE first,testcpu
	200
	201	REAL SSUM
	202	REAL temps0,temps1,temps2,temps3,temps4,temps5
	203	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	204
	205
	206	DATA first,testcpu/.true.,.false./
	207
	208	IF(first) THEN
	209	temps1=0.
	210	temps2=0.
	211	temps3=0.
	212	temps4=0.
	213	temps5=0.
	214	first=.false.
	215	ENDIF
	216
	217	c calcul de la pente a droite et a gauche de la maille
	218
	219
	220	IF (pente_max.gt.-1.e-5) THEN
	221	c IF (pente_max.gt.10) THEN
	222
	223	c calcul des pentes avec limitation, Van Leer scheme I:
	224	c -----------------------------------------------------
	225
	226	c calcul de la pente aux points u
	227	DO l = 1, llm
	228	DO ij=iip2,ip1jm-1
	229	dxqu(ij)=q(ij+1,l)-q(ij,l)
	230	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	231	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	232	ENDDO
	233	DO ij=iip1+iip1,ip1jm,iip1
	234	dxqu(ij)=dxqu(ij-iim)
	235	c sigu(ij)=sigu(ij-iim)
	236	ENDDO
	237
	238	DO ij=iip2,ip1jm
	239	adxqu(ij)=abs(dxqu(ij))
	240	ENDDO
	241
	242	c calcul de la pente maximum dans la maille en valeur absolue
	243
	244	DO ij=iip2+1,ip1jm
	245	dxqmax(ij,l)=pente_max*
	246	, min(adxqu(ij-1),adxqu(ij))
	247	c limitation subtile
	248	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	249
	250
	251	ENDDO
	252
	253	DO ij=iip1+iip1,ip1jm,iip1
	254	dxqmax(ij-iim,l)=dxqmax(ij,l)
	255	ENDDO
	256
	257	DO ij=iip2+1,ip1jm
	258	#ifdef CRAY
	259	dxq(ij,l)=
	260	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	261	#else
	262	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	263	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	264	ELSE
	265	c extremum local
	266	dxq(ij,l)=0.
	267	ENDIF
	268	#endif
	269	dxq(ij,l)=0.5*dxq(ij,l)
	270	dxq(ij,l)=
	271	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	272	ENDDO
	273
	274	ENDDO ! l=1,llm
	275
	276	ELSE ! (pente_max.lt.-1.e-5)
	277
	278	c Pentes produits:
	279	c ----------------
	280
	281	DO l = 1, llm
	282	DO ij=iip2,ip1jm-1
	283	dxqu(ij)=q(ij+1,l)-q(ij,l)
	284	ENDDO
	285	DO ij=iip1+iip1,ip1jm,iip1
	286	dxqu(ij)=dxqu(ij-iim)
	287	ENDDO
	288
	289	DO ij=iip2+1,ip1jm
	290	zz(ij)=dxqu(ij-1)*dxqu(ij)
	291	zz(ij)=zz(ij)+zz(ij)
	292	IF(zz(ij).gt.0) THEN
	293	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	294	ELSE
	295	c extremum local
	296	dxq(ij,l)=0.
	297	ENDIF
	298	ENDDO
	299
	300	ENDDO
	301
	302	ENDIF ! (pente_max.lt.-1.e-5)
	303
	304	c bouclage de la pente en iip1:
	305	c -----------------------------
	306
	307	DO l=1,llm
	308	DO ij=iip1+iip1,ip1jm,iip1
	309	dxq(ij-iim,l)=dxq(ij,l)
	310	ENDDO
	311
	312	DO ij=1,ip1jmp1
	313	iadvplus(ij,l)=0
	314	ENDDO
	315
	316	ENDDO
	317
	318
	319	c calcul des flux a gauche et a droite
	320
	321	#ifdef CRAY
	322	c--pas encore modification sur Qsat
	323	DO l=1,llm
	324	DO ij=iip2,ip1jm-1
	325	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	326	, 1.+u_m(ij,l)/masse(ij+1,l),
	327	, u_m(ij,l))
	328	zdum(ij,l)=0.5*zdum(ij,l)
	329	u_mq(ij,l)=cvmgp(
	330	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	331	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	332	, u_m(ij,l))
	333	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	334	ENDDO
	335	ENDDO
	336	#else
	337	c on cumule le flux correspondant a toutes les mailles dont la masse
	338	c au travers de la paroi pENDant le pas de temps.
	339	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
	340	DO l=1,llm
	341	DO ij=iip2,ip1jm-1
	342	IF (u_m(ij,l).gt.0.) THEN
	343	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	344	u_mq(ij,l)=u_m(ij,l)*
	345	$ min(q(ij,l)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
	346	ELSE
	347	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	348	u_mq(ij,l)=u_m(ij,l)*
	349	$ min(q(ij+1,l)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
	350	ENDIF
	351	ENDDO
	352	ENDDO
	353	#endif
	354
	355
	356	c detection des points ou on advecte plus que la masse de la
	357	c maille
	358	DO l=1,llm
	359	DO ij=iip2,ip1jm-1
	360	IF(zdum(ij,l).lt.0) THEN
	361	iadvplus(ij,l)=1
	362	u_mq(ij,l)=0.
	363	ENDIF
	364	ENDDO
	365	ENDDO
	366	DO l=1,llm
	367	DO ij=iip1+iip1,ip1jm,iip1
	368	iadvplus(ij,l)=iadvplus(ij-iim,l)
	369	ENDDO
	370	ENDDO
	371
	372
	373
	374	c traitement special pour le cas ou on advecte en longitude plus que le
	375	c contenu de la maille.
	376	c cette partie est mal vectorisee.
	377
	378	c pas d'influence de la pression saturante (pour l'instant)
	379
	380	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	381
	382	n0=0
	383	DO l=1,llm
	384	nl(l)=0
	385	DO ij=iip2,ip1jm
	386	nl(l)=nl(l)+iadvplus(ij,l)
	387	ENDDO
	388	n0=n0+nl(l)
	389	ENDDO
	390
[595]	391	IF(n0.gt.0) THEN
[524]	392	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
	393	ccc & ,'contenu de la maille : ',n0
	394
	395	DO l=1,llm
	396	IF(nl(l).gt.0) THEN
	397	iju=0
	398	c indicage des mailles concernees par le traitement special
	399	DO ij=iip2,ip1jm
	400	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	401	iju=iju+1
	402	indu(iju)=ij
	403	ENDIF
	404	ENDDO
	405	niju=iju
	406	c PRINT*,'niju,nl',niju,nl(l)
	407
	408	c traitement des mailles
	409	DO iju=1,niju
	410	ij=indu(iju)
	411	j=(ij-1)/iip1+1
	412	zu_m=u_m(ij,l)
	413	u_mq(ij,l)=0.
	414	IF(zu_m.gt.0.) THEN
	415	ijq=ij
	416	i=ijq-(j-1)*iip1
	417	c accumulation pour les mailles completements advectees
	418	do while(zu_m.gt.masse(ijq,l))
	419	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	420	zu_m=zu_m-masse(ijq,l)
	421	i=mod(i-2+iim,iim)+1
	422	ijq=(j-1)*iip1+i
	423	ENDDO
	424	c ajout de la maille non completement advectee
	425	u_mq(ij,l)=u_mq(ij,l)+zu_m*
	426	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
	427	ELSE
	428	ijq=ij+1
	429	i=ijq-(j-1)*iip1
	430	c accumulation pour les mailles completements advectees
	431	do while(-zu_m.gt.masse(ijq,l))
	432	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	433	zu_m=zu_m+masse(ijq,l)
	434	i=mod(i,iim)+1
	435	ijq=(j-1)*iip1+i
	436	ENDDO
	437	c ajout de la maille non completement advectee
	438	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
	439	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
	440	ENDIF
	441	ENDDO
	442	ENDIF
	443	ENDDO
	444	ENDIF ! n0.gt.0
	445
	446
	447
	448	c bouclage en latitude
	449
	450	DO l=1,llm
	451	DO ij=iip1+iip1,ip1jm,iip1
	452	u_mq(ij,l)=u_mq(ij-iim,l)
	453	ENDDO
	454	ENDDO
	455
	456
	457	c calcul des tendances
	458
	459	DO l=1,llm
	460	DO ij=iip2+1,ip1jm
	461	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	462	q(ij,l)=(q(ij,l)*masse(ij,l)+
	463	& u_mq(ij-1,l)-u_mq(ij,l))
	464	& /new_m
	465	masse(ij,l)=new_m
	466	ENDDO
	467	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	468	DO ij=iip1+iip1,ip1jm,iip1
	469	q(ij-iim,l)=q(ij,l)
	470	masse(ij-iim,l)=masse(ij,l)
	471	ENDDO
	472	ENDDO
	473
	474	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	475	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	476
	477
	478	RETURN
	479	END
	480	SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
	481	c
	482	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	483	c
	484	c ********************************************************************
	485	c Shema d'advection " pseudo amont " .
	486	c ********************************************************************
	487	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	488	c qsat est un argument de sortie pour le s-pg ....
	489	c
	490	c
	491	c --------------------------------------------------------------------
	492	IMPLICIT NONE
	493	c
	494	#include "dimensions.h"
	495	#include "paramet.h"
	496	#include "logic.h"
	497	#include "comvert.h"
	498	#include "comconst.h"
	499	#include "comgeom.h"
	500	c
	501	c
	502	c Arguments:
	503	c ----------
	504	REAL masse(ip1jmp1,llm),pente_max
	505	REAL masse_adv_v( ip1jm,llm)
	506	REAL q(ip1jmp1,llm)
	507	REAL qsat(ip1jmp1,llm)
	508	c
	509	c Local
	510	c ---------
	511	c
	512	INTEGER i,ij,l
	513	c
	514	REAL airej2,airejjm,airescb(iim),airesch(iim)
	515	REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
	516	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
	517	REAL qbyv(ip1jm,llm)
	518
	519	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	520	c REAL newq,oldmasse
	521	Logical first,testcpu
	522	REAL temps0,temps1,temps2,temps3,temps4,temps5
	523	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	524	SAVE first,testcpu
	525
	526	REAL convpn,convps,convmpn,convmps
	527	REAL sinlon(iip1),sinlondlon(iip1)
	528	REAL coslon(iip1),coslondlon(iip1)
	529	SAVE sinlon,coslon,sinlondlon,coslondlon
	530	SAVE airej2,airejjm
	531	c
	532	c
	533	REAL SSUM
	534
	535	DATA first,testcpu/.true.,.false./
	536	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
	537
	538	IF(first) THEN
	539	PRINT*,'Shema Amont nouveau appele dans Vanleer '
	540	first=.false.
	541	do i=2,iip1
	542	coslon(i)=cos(rlonv(i))
	543	sinlon(i)=sin(rlonv(i))
	544	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	545	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	546	ENDDO
	547	coslon(1)=coslon(iip1)
	548	coslondlon(1)=coslondlon(iip1)
	549	sinlon(1)=sinlon(iip1)
	550	sinlondlon(1)=sinlondlon(iip1)
	551	airej2 = SSUM( iim, aire(iip2), 1 )
	552	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
	553	ENDIF
	554
	555	c
	556
	557
	558	DO l = 1, llm
	559	c
	560	c --------------------------------
	561	c CALCUL EN LATITUDE
	562	c --------------------------------
	563
	564	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
	565	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
	566	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
	567
	568	DO i = 1, iim
	569	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
	570	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
	571	ENDDO
	572	qpns = SSUM( iim, airescb ,1 ) / airej2
	573	qpsn = SSUM( iim, airesch ,1 ) / airejjm
	574
	575	c calcul des pentes aux points v
	576
	577	DO ij=1,ip1jm
	578	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
	579	adyqv(ij)=abs(dyqv(ij))
	580	ENDDO
	581
	582	c calcul des pentes aux points scalaires
	583
	584	DO ij=iip2,ip1jm
	585	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
	586	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
	587	dyqmax(ij)=pente_max*dyqmax(ij)
	588	ENDDO
	589
	590	c calcul des pentes aux poles
	591
	592	DO ij=1,iip1
	593	dyq(ij,l)=qpns-q(ij+iip1,l)
	594	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
	595	ENDDO
	596
	597	c filtrage de la derivee
	598	dyn1=0.
	599	dys1=0.
	600	dyn2=0.
	601	dys2=0.
	602	DO ij=1,iim
	603	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
	604	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
	605	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
	606	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
	607	ENDDO
	608	DO ij=1,iip1
	609	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
	610	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
	611	ENDDO
	612
	613	c calcul des pentes limites aux poles
	614
	615	fn=1.
	616	fs=1.
	617	DO ij=1,iim
	618	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
	619	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
	620	ENDIF
	621	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
	622	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
	623	ENDIF
	624	ENDDO
	625	DO ij=1,iip1
	626	dyq(ij,l)=fn*dyq(ij,l)
	627	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
	628	ENDDO
	629
	630	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	631	C En memoire de dIFferents tests sur la
	632	C limitation des pentes aux poles.
	633	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	634	C PRINT*,dyq(1)
	635	C PRINT*,dyqv(iip1+1)
	636	C apn=abs(dyq(1)/dyqv(iip1+1))
	637	C PRINT*,dyq(ip1jm+1)
	638	C PRINT*,dyqv(ip1jm-iip1+1)
	639	C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
	640	C DO ij=2,iim
	641	C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
	642	C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
	643	C ENDDO
	644	C apn=min(pente_max/apn,1.)
	645	C aps=min(pente_max/aps,1.)
	646	C
	647	C
	648	C cas ou on a un extremum au pole
	649	C
	650	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	651	C & apn=0.
	652	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	653	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	654	C & aps=0.
	655	C
	656	C limitation des pentes aux poles
	657	C DO ij=1,iip1
	658	C dyq(ij)=apn*dyq(ij)
	659	C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
	660	C ENDDO
	661	C
	662	C test
	663	C DO ij=1,iip1
	664	C dyq(iip1+ij)=0.
	665	C dyq(ip1jm+ij-iip1)=0.
	666	C ENDDO
	667	C DO ij=1,ip1jmp1
	668	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
	669	C ENDDO
	670	C
	671	C changement 10 07 96
	672	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	673	C & THEN
	674	C DO ij=1,iip1
	675	C dyqmax(ij)=0.
	676	C ENDDO
	677	C ELSE
	678	C DO ij=1,iip1
	679	C dyqmax(ij)=pente_max*abs(dyqv(ij))
	680	C ENDDO
	681	C ENDIF
	682	C
	683	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	684	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	685	C &THEN
	686	C DO ij=ip1jm+1,ip1jmp1
	687	C dyqmax(ij)=0.
	688	C ENDDO
	689	C ELSE
	690	C DO ij=ip1jm+1,ip1jmp1
	691	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
	692	C ENDDO
	693	C ENDIF
	694	C fin changement 10 07 96
	695	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	696
	697	c calcul des pentes limitees
	698
	699	DO ij=iip2,ip1jm
	700	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
	701	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
	702	ELSE
	703	dyq(ij,l)=0.
	704	ENDIF
	705	ENDDO
	706
	707	ENDDO
	708
	709	DO l=1,llm
	710	DO ij=1,ip1jm
	711	IF( masse_adv_v(ij,l).GT.0. ) THEN
	712	qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) +
	713	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
	714	ELSE
	715	qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
	716	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
	717	ENDIF
	718	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
	719	ENDDO
	720	ENDDO
	721
	722
	723	DO l=1,llm
	724	DO ij=iip2,ip1jm
	725	newmasse=masse(ij,l)
	726	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
	727	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
	728	& /newmasse
	729	masse(ij,l)=newmasse
	730	ENDDO
	731	c.-. ancienne version
	732	convpn=SSUM(iim,qbyv(1,l),1)/apoln
	733	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
	734	DO ij = 1,iip1
	735	newmasse=masse(ij,l)+convmpn*aire(ij)
	736	q(ij,l)=(q(ij,l)masse(ij,l)+convpnaire(ij))/
	737	& newmasse
	738	masse(ij,l)=newmasse
	739	ENDDO
	740	convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
	741	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
	742	DO ij = ip1jm+1,ip1jmp1
	743	newmasse=masse(ij,l)+convmps*aire(ij)
	744	q(ij,l)=(q(ij,l)masse(ij,l)+convpsaire(ij))/
	745	& newmasse
	746	masse(ij,l)=newmasse
	747	ENDDO
	748	c.-. fin ancienne version
	749
	750	c._. nouvelle version
	751	c convpn=SSUM(iim,qbyv(1,l),1)
	752	c convmpn=ssum(iim,masse_adv_v(1,l),1)
	753	c oldmasse=ssum(iim,masse(1,l),1)
	754	c newmasse=oldmasse+convmpn
	755	c newq=(q(1,l)*oldmasse+convpn)/newmasse
	756	c newmasse=newmasse/apoln
	757	c DO ij = 1,iip1
	758	c q(ij,l)=newq
	759	c masse(ij,l)=newmasse*aire(ij)
	760	c ENDDO
	761	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	762	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	763	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
	764	c newmasse=oldmasse+convmps
	765	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
	766	c newmasse=newmasse/apols
	767	c DO ij = ip1jm+1,ip1jmp1
	768	c q(ij,l)=newq
	769	c masse(ij,l)=newmasse*aire(ij)
	770	c ENDDO
	771	c._. fin nouvelle version
	772	ENDDO
	773
	774	RETURN
	775	END

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

Download in other formats:

Original Format