Context Navigation

vlspltqs.F @ 5428

Last change on this file since 5428 was 4482, checked in by lguez, 21 months ago
Sync latest trunk changes to branch LMDZ_ECRad
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 24.4 KB

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

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

Download in other formats:

Original Format