Context Navigation

vlspltqs.F @ 1243

Last change on this file since 1243 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 21.1 KB

Rev	Line
[135]	1	SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
	2	, p,pk,teta )
	3	c
	4	c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
	5	c
	6	c ********************************************************************
	7	c Shema d'advection " pseudo amont " .
	8	c + test sur humidite specifique: Q advecte< Qsat aval
	9	c (F. Codron, 10/99)
	10	c ********************************************************************
	11	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	12	c
	13	c pente_max facteur de limitation des pentes: 2 en general
	14	c 0 pour un schema amont
	15	c pbaru,pbarv,w flux de masse en u ,v ,w
	16	c pdt pas de temps
	17	c
	18	c teta temperature potentielle, p pression aux interfaces,
	19	c pk exner au milieu des couches necessaire pour calculer Qsat
	20	c --------------------------------------------------------------------
	21	IMPLICIT NONE
	22	c
	23	#include "dimensions.h"
	24	#include "paramet.h"
	25	#include "logic.h"
	26	#include "comvert.h"
	27	#include "comconst.h"
	28
	29	c
	30	c Arguments:
	31	c ----------
	32	REAL masse(ip1jmp1,llm),pente_max
	33	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
	34	REAL q(ip1jmp1,llm)
	35	REAL w(ip1jmp1,llm),pdt
	36	REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
	37	c
	38	c Local
	39	c ---------
	40	c
	41	INTEGER i,ij,l,j,ii
	42	c
	43	REAL qsat(ip1jmp1,llm)
	44	REAL zm(ip1jmp1,llm)
	45	REAL mu(ip1jmp1,llm)
	46	REAL mv(ip1jm,llm)
	47	REAL mw(ip1jmp1,llm+1)
	48	REAL zq(ip1jmp1,llm)
	49	REAL temps1,temps2,temps3
	50	REAL zzpbar, zzw
	51	LOGICAL testcpu
	52	SAVE testcpu
	53	SAVE temps1,temps2,temps3
	54
	55	REAL qmin,qmax
	56	DATA qmin,qmax/0.,1.e33/
	57	DATA testcpu/.false./
	58	DATA temps1,temps2,temps3/0.,0.,0./
	59
	60	c--pour rapport de melange saturant--
	61
	62	REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
	63	REAL ptarg,pdelarg,foeew,zdelta
	64	REAL tempe(ip1jmp1)
	65
	66	c fonction psat(T)
	67
	68	FOEEW ( PTARG,PDELARG ) = EXP (
	69	* (R3LES(1.-PDELARG)+R3IESPDELARG) * (PTARG-RTT)
	70	* / (PTARG-(R4LES(1.-PDELARG)+R4IESPDELARG)) )
	71
	72	r2es = 380.11733
	73	r3les = 17.269
	74	r3ies = 21.875
	75	r4les = 35.86
	76	r4ies = 7.66
	77	retv = 0.6077667
	78	rtt = 273.16
	79
	80	c-- Calcul de Qsat en chaque point
	81	c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
	82	c pour eviter une exponentielle.
	83	DO l = 1, llm
	84	DO ij = 1, ip1jmp1
	85	tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
	86	ENDDO
	87	DO ij = 1, ip1jmp1
	88	zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
	89	play = 0.5*(p(ij,l)+p(ij,l+1))
	90	qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
	91	qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
	92	ENDDO
	93	ENDDO
	94
	95
	96	c PRINT*,'Debut vlsplt version debug sans vlyqs'
	97
	98	zzpbar = 0.5 * pdt
	99	zzw = pdt
	100	DO l=1,llm
	101	DO ij = iip2,ip1jm
	102	mu(ij,l)=pbaru(ij,l) * zzpbar
	103	ENDDO
	104	DO ij=1,ip1jm
	105	mv(ij,l)=pbarv(ij,l) * zzpbar
	106	ENDDO
	107	DO ij=1,ip1jmp1
	108	mw(ij,l)=w(ij,l) * zzw
	109	ENDDO
	110	ENDDO
	111
	112	DO ij=1,ip1jmp1
	113	mw(ij,llm+1)=0.
	114	ENDDO
	115
	116	CALL SCOPY(ijp1llm,q,1,zq,1)
	117	CALL SCOPY(ijp1llm,masse,1,zm,1)
	118
	119	c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
	120	call vlxqs(zq,pente_max,zm,mu,qsat)
	121
	122
	123	c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
	124
	125	call vlyqs(zq,pente_max,zm,mv,qsat)
	126
	127
	128	c call minmaxq(zq,qmin,qmax,'avant vlz ')
	129
	130	call vlz(zq,pente_max,zm,mw)
	131
	132
	133	c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
	134	c call minmaxq(zm,qmin,qmax,'M avant vlyqs ')
	135
	136	call vlyqs(zq,pente_max,zm,mv,qsat)
	137
	138
	139	c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
	140	c call minmaxq(zm,qmin,qmax,'M avant vlxqs ')
	141
	142	call vlxqs(zq,pente_max,zm,mu,qsat)
	143
	144	c call minmaxq(zq,qmin,qmax,'apres vlxqs ')
	145	c call minmaxq(zm,qmin,qmax,'M apres vlxqs ')
	146
	147
	148	DO l=1,llm
	149	DO ij=1,ip1jmp1
	150	q(ij,l)=zq(ij,l)
	151	ENDDO
	152	DO ij=1,ip1jm+1,iip1
	153	q(ij+iim,l)=q(ij,l)
	154	ENDDO
	155	ENDDO
	156
	157	RETURN
	158	END
	159	SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
	160	c
	161	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	162	c
	163	c ********************************************************************
	164	c Shema d'advection " pseudo amont " .
	165	c ********************************************************************
	166	c
	167	c --------------------------------------------------------------------
	168	IMPLICIT NONE
	169	c
	170	#include "dimensions.h"
	171	#include "paramet.h"
	172	#include "logic.h"
	173	#include "comvert.h"
	174	#include "comconst.h"
	175	c
	176	c
	177	c Arguments:
	178	c ----------
	179	REAL masse(ip1jmp1,llm),pente_max
	180	REAL u_m( ip1jmp1,llm )
	181	REAL q(ip1jmp1,llm)
	182	REAL qsat(ip1jmp1,llm)
	183	c
	184	c Local
	185	c ---------
	186	c
	187	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
	188	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
	189	c
	190	REAL new_m,zu_m,zdum(ip1jmp1,llm)
	191	REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
	192	REAL zz(ip1jmp1)
	193	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
	194	REAL u_mq(ip1jmp1,llm)
	195
	196	Logical first,testcpu
	197	SAVE first,testcpu
	198
	199	REAL SSUM
	200	EXTERNAL SSUM
	201	REAL temps0,temps1,temps2,temps3,temps4,temps5
	202	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	203
	204
	205	DATA first,testcpu/.true.,.false./
	206
	207	IF(first) THEN
	208	temps1=0.
	209	temps2=0.
	210	temps3=0.
	211	temps4=0.
	212	temps5=0.
	213	first=.false.
	214	ENDIF
	215
	216	c calcul de la pente a droite et a gauche de la maille
	217
	218
	219	IF (pente_max.gt.-1.e-5) THEN
	220	c IF (pente_max.gt.10) THEN
	221
	222	c calcul des pentes avec limitation, Van Leer scheme I:
	223	c -----------------------------------------------------
	224
	225	c calcul de la pente aux points u
	226	DO l = 1, llm
	227	DO ij=iip2,ip1jm-1
	228	dxqu(ij)=q(ij+1,l)-q(ij,l)
	229	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	230	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	231	ENDDO
	232	DO ij=iip1+iip1,ip1jm,iip1
	233	dxqu(ij)=dxqu(ij-iim)
	234	c sigu(ij)=sigu(ij-iim)
	235	ENDDO
	236
	237	DO ij=iip2,ip1jm
	238	adxqu(ij)=abs(dxqu(ij))
	239	ENDDO
	240
	241	c calcul de la pente maximum dans la maille en valeur absolue
	242
	243	DO ij=iip2+1,ip1jm
	244	dxqmax(ij,l)=pente_max*
	245	, min(adxqu(ij-1),adxqu(ij))
	246	c limitation subtile
	247	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	248
	249
	250	ENDDO
	251
	252	DO ij=iip1+iip1,ip1jm,iip1
	253	dxqmax(ij-iim,l)=dxqmax(ij,l)
	254	ENDDO
	255
	256	DO ij=iip2+1,ip1jm
	257	#ifdef CRAY
	258	dxq(ij,l)=
	259	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	260	#else
	261	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	262	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	263	ELSE
	264	c extremum local
	265	dxq(ij,l)=0.
	266	ENDIF
	267	#endif
	268	dxq(ij,l)=0.5*dxq(ij,l)
	269	dxq(ij,l)=
	270	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	271	ENDDO
	272
	273	ENDDO ! l=1,llm
	274
	275	ELSE ! (pente_max.lt.-1.e-5)
	276
	277	c Pentes produits:
	278	c ----------------
	279
	280	DO l = 1, llm
	281	DO ij=iip2,ip1jm-1
	282	dxqu(ij)=q(ij+1,l)-q(ij,l)
	283	ENDDO
	284	DO ij=iip1+iip1,ip1jm,iip1
	285	dxqu(ij)=dxqu(ij-iim)
	286	ENDDO
	287
	288	DO ij=iip2+1,ip1jm
	289	zz(ij)=dxqu(ij-1)*dxqu(ij)
	290	zz(ij)=zz(ij)+zz(ij)
	291	IF(zz(ij).gt.0) THEN
	292	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	293	ELSE
	294	c extremum local
	295	dxq(ij,l)=0.
	296	ENDIF
	297	ENDDO
	298
	299	ENDDO
	300
	301	ENDIF ! (pente_max.lt.-1.e-5)
	302
	303	c bouclage de la pente en iip1:
	304	c -----------------------------
	305
	306	DO l=1,llm
	307	DO ij=iip1+iip1,ip1jm,iip1
	308	dxq(ij-iim,l)=dxq(ij,l)
	309	ENDDO
	310
	311	DO ij=1,ip1jmp1
	312	iadvplus(ij,l)=0
	313	ENDDO
	314
	315	ENDDO
	316
	317
	318	c calcul des flux a gauche et a droite
	319
	320	#ifdef CRAY
	321	c--pas encore modification sur Qsat
	322	DO l=1,llm
	323	DO ij=iip2,ip1jm-1
	324	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	325	, 1.+u_m(ij,l)/masse(ij+1,l),
	326	, u_m(ij,l))
	327	zdum(ij,l)=0.5*zdum(ij,l)
	328	u_mq(ij,l)=cvmgp(
	329	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	330	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	331	, u_m(ij,l))
	332	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	333	ENDDO
	334	ENDDO
	335	#else
	336	c on cumule le flux correspondant a toutes les mailles dont la masse
	337	c au travers de la paroi pENDant le pas de temps.
	338	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
	339	DO l=1,llm
	340	DO ij=iip2,ip1jm-1
	341	IF (u_m(ij,l).gt.0.) THEN
	342	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	343	u_mq(ij,l)=u_m(ij,l)*
	344	$ min(q(ij,l)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
	345	ELSE
	346	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	347	u_mq(ij,l)=u_m(ij,l)*
	348	$ min(q(ij+1,l)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
	349	ENDIF
	350	ENDDO
	351	ENDDO
	352	#endif
	353
	354
	355	c detection des points ou on advecte plus que la masse de la
	356	c maille
	357	DO l=1,llm
	358	DO ij=iip2,ip1jm-1
	359	IF(zdum(ij,l).lt.0) THEN
	360	iadvplus(ij,l)=1
	361	u_mq(ij,l)=0.
	362	ENDIF
	363	ENDDO
	364	ENDDO
	365	DO l=1,llm
	366	DO ij=iip1+iip1,ip1jm,iip1
	367	iadvplus(ij,l)=iadvplus(ij-iim,l)
	368	ENDDO
	369	ENDDO
	370
	371
	372
	373	c traitement special pour le cas ou on advecte en longitude plus que le
	374	c contenu de la maille.
	375	c cette partie est mal vectorisee.
	376
	377	c pas d'influence de la pression saturante (pour l'instant)
	378
	379	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	380
	381	n0=0
	382	DO l=1,llm
	383	nl(l)=0
	384	DO ij=iip2,ip1jm
	385	nl(l)=nl(l)+iadvplus(ij,l)
	386	ENDDO
	387	n0=n0+nl(l)
	388	ENDDO
	389
	390	IF(n0.gt.1) THEN
	391	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
	392	ccc & ,'contenu de la maille : ',n0
	393
	394	DO l=1,llm
	395	IF(nl(l).gt.0) THEN
	396	iju=0
	397	c indicage des mailles concernees par le traitement special
	398	DO ij=iip2,ip1jm
	399	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	400	iju=iju+1
	401	indu(iju)=ij
	402	ENDIF
	403	ENDDO
	404	niju=iju
	405	c PRINT*,'niju,nl',niju,nl(l)
	406
	407	c traitement des mailles
	408	DO iju=1,niju
	409	ij=indu(iju)
	410	j=(ij-1)/iip1+1
	411	zu_m=u_m(ij,l)
	412	u_mq(ij,l)=0.
	413	IF(zu_m.gt.0.) THEN
	414	ijq=ij
	415	i=ijq-(j-1)*iip1
	416	c accumulation pour les mailles completements advectees
	417	do while(zu_m.gt.masse(ijq,l))
	418	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	419	zu_m=zu_m-masse(ijq,l)
	420	i=mod(i-2+iim,iim)+1
	421	ijq=(j-1)*iip1+i
	422	ENDDO
	423	c ajout de la maille non completement advectee
	424	u_mq(ij,l)=u_mq(ij,l)+zu_m*
	425	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
	426	ELSE
	427	ijq=ij+1
	428	i=ijq-(j-1)*iip1
	429	c accumulation pour les mailles completements advectees
	430	do while(-zu_m.gt.masse(ijq,l))
	431	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	432	zu_m=zu_m+masse(ijq,l)
	433	i=mod(i,iim)+1
	434	ijq=(j-1)*iip1+i
	435	ENDDO
	436	c ajout de la maille non completement advectee
	437	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
	438	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
	439	ENDIF
	440	ENDDO
	441	ENDIF
	442	ENDDO
	443	ENDIF ! n0.gt.0
	444
	445
	446
	447	c bouclage en latitude
	448
	449	DO l=1,llm
	450	DO ij=iip1+iip1,ip1jm,iip1
	451	u_mq(ij,l)=u_mq(ij-iim,l)
	452	ENDDO
	453	ENDDO
	454
	455
	456	c calcul des tendances
	457
	458	DO l=1,llm
	459	DO ij=iip2+1,ip1jm
	460	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	461	q(ij,l)=(q(ij,l)*masse(ij,l)+
	462	& u_mq(ij-1,l)-u_mq(ij,l))
	463	& /new_m
	464	masse(ij,l)=new_m
	465	ENDDO
	466	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	467	DO ij=iip1+iip1,ip1jm,iip1
	468	q(ij-iim,l)=q(ij,l)
	469	masse(ij-iim,l)=masse(ij,l)
	470	ENDDO
	471	ENDDO
	472
	473	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	474	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	475
	476
	477	RETURN
	478	END
	479	SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
	480	c
	481	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	482	c
	483	c ********************************************************************
	484	c Shema d'advection " pseudo amont " .
	485	c ********************************************************************
	486	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	487	c qsat est un argument de sortie pour le s-pg ....
	488	c
	489	c
	490	c --------------------------------------------------------------------
	491	IMPLICIT NONE
	492	c
	493	#include "dimensions.h"
	494	#include "paramet.h"
	495	#include "logic.h"
	496	#include "comvert.h"
	497	#include "comconst.h"
	498	#include "comgeom.h"
	499	c
	500	c
	501	c Arguments:
	502	c ----------
	503	REAL masse(ip1jmp1,llm),pente_max
	504	REAL masse_adv_v( ip1jm,llm)
	505	REAL q(ip1jmp1,llm)
	506	REAL qsat(ip1jmp1,llm)
	507	c
	508	c Local
	509	c ---------
	510	c
	511	INTEGER i,ij,l
	512	c
	513	REAL airej2,airejjm,airescb(iim),airesch(iim)
	514	REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
	515	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
	516	REAL qbyv(ip1jm,llm)
	517
	518	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	519	c REAL newq,oldmasse
	520	Logical first,testcpu
	521	REAL temps0,temps1,temps2,temps3,temps4,temps5
	522	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	523	SAVE first,testcpu
	524
	525	REAL convpn,convps,convmpn,convmps
	526	REAL sinlon(iip1),sinlondlon(iip1)
	527	REAL coslon(iip1),coslondlon(iip1)
	528	SAVE sinlon,coslon,sinlondlon,coslondlon
	529	SAVE airej2,airejjm
	530	c
	531	c
	532	REAL SSUM
	533	EXTERNAL 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