Context Navigation

vlspltqs.F @ 5453

Last change on this file since 5453 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
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.6 KB

Rev	Line
[1520]	1	c
	2	c $Id: vlspltqs.F 4368 2022-12-05 23:01:16Z aborella $
	3	c
[524]	4	SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
[2270]	5	, p,pk,teta,iq )
[4368]	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)
[4368]	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)
[4368]	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
[4368]	162	do ifils=1,tracers(iq)%nqDescen
	163	iq2=tracers(iq)%iqDescen(ifils)
[2270]	164	DO l=1,llm
[4368]	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)
[4368]	170	ENDDO
[2270]	171	ENDDO
[4368]	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)
[4368]	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	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
[2270]	251	c sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
[524]	252	ENDDO
	253	DO ij=iip1+iip1,ip1jm,iip1
	254	dxqu(ij)=dxqu(ij-iim)
	255	c sigu(ij)=sigu(ij-iim)
	256	ENDDO
	257
	258	DO ij=iip2,ip1jm
	259	adxqu(ij)=abs(dxqu(ij))
	260	ENDDO
	261
	262	c calcul de la pente maximum dans la maille en valeur absolue
	263
	264	DO ij=iip2+1,ip1jm
	265	dxqmax(ij,l)=pente_max*
	266	, min(adxqu(ij-1),adxqu(ij))
	267	c limitation subtile
	268	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	269
	270
	271	ENDDO
	272
	273	DO ij=iip1+iip1,ip1jm,iip1
	274	dxqmax(ij-iim,l)=dxqmax(ij,l)
	275	ENDDO
	276
	277	DO ij=iip2+1,ip1jm
	278	#ifdef CRAY
	279	dxq(ij,l)=
	280	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	281	#else
	282	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	283	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	284	ELSE
	285	c extremum local
	286	dxq(ij,l)=0.
	287	ENDIF
	288	#endif
	289	dxq(ij,l)=0.5*dxq(ij,l)
	290	dxq(ij,l)=
	291	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	292	ENDDO
	293
	294	ENDDO ! l=1,llm
	295
	296	ELSE ! (pente_max.lt.-1.e-5)
	297
	298	c Pentes produits:
	299	c ----------------
	300
	301	DO l = 1, llm
	302	DO ij=iip2,ip1jm-1
[2270]	303	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
[524]	304	ENDDO
	305	DO ij=iip1+iip1,ip1jm,iip1
	306	dxqu(ij)=dxqu(ij-iim)
	307	ENDDO
	308
	309	DO ij=iip2+1,ip1jm
	310	zz(ij)=dxqu(ij-1)*dxqu(ij)
	311	zz(ij)=zz(ij)+zz(ij)
	312	IF(zz(ij).gt.0) THEN
	313	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	314	ELSE
	315	c extremum local
	316	dxq(ij,l)=0.
	317	ENDIF
	318	ENDDO
	319
	320	ENDDO
	321
	322	ENDIF ! (pente_max.lt.-1.e-5)
	323
	324	c bouclage de la pente en iip1:
	325	c -----------------------------
	326
	327	DO l=1,llm
	328	DO ij=iip1+iip1,ip1jm,iip1
	329	dxq(ij-iim,l)=dxq(ij,l)
	330	ENDDO
	331
	332	DO ij=1,ip1jmp1
	333	iadvplus(ij,l)=0
	334	ENDDO
	335
	336	ENDDO
	337
	338
	339	c calcul des flux a gauche et a droite
	340
	341	#ifdef CRAY
	342	c--pas encore modification sur Qsat
	343	DO l=1,llm
	344	DO ij=iip2,ip1jm-1
[2270]	345	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
	346	, 1.+u_m(ij,l)/masse(ij+1,l,iq),
[524]	347	, u_m(ij,l))
	348	zdum(ij,l)=0.5*zdum(ij,l)
	349	u_mq(ij,l)=cvmgp(
[2270]	350	, q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
	351	, q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
[524]	352	, u_m(ij,l))
	353	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	354	ENDDO
	355	ENDDO
	356	#else
	357	c on cumule le flux correspondant a toutes les mailles dont la masse
	358	c au travers de la paroi pENDant le pas de temps.
	359	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
	360	DO l=1,llm
	361	DO ij=iip2,ip1jm-1
	362	IF (u_m(ij,l).gt.0.) THEN
[2270]	363	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
[524]	364	u_mq(ij,l)=u_m(ij,l)*
[2270]	365	$ min(q(ij,l,iq)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
[524]	366	ELSE
[2270]	367	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
[524]	368	u_mq(ij,l)=u_m(ij,l)*
[2270]	369	$ min(q(ij+1,l,iq)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
[524]	370	ENDIF
	371	ENDDO
	372	ENDDO
	373	#endif
	374
	375
	376	c detection des points ou on advecte plus que la masse de la
	377	c maille
	378	DO l=1,llm
	379	DO ij=iip2,ip1jm-1
	380	IF(zdum(ij,l).lt.0) THEN
	381	iadvplus(ij,l)=1
	382	u_mq(ij,l)=0.
	383	ENDIF
	384	ENDDO
	385	ENDDO
	386	DO l=1,llm
	387	DO ij=iip1+iip1,ip1jm,iip1
	388	iadvplus(ij,l)=iadvplus(ij-iim,l)
	389	ENDDO
	390	ENDDO
	391
	392
	393
	394	c traitement special pour le cas ou on advecte en longitude plus que le
	395	c contenu de la maille.
	396	c cette partie est mal vectorisee.
	397
	398	c pas d'influence de la pression saturante (pour l'instant)
	399
	400	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	401
	402	n0=0
	403	DO l=1,llm
	404	nl(l)=0
	405	DO ij=iip2,ip1jm
	406	nl(l)=nl(l)+iadvplus(ij,l)
	407	ENDDO
	408	n0=n0+nl(l)
	409	ENDDO
	410
[595]	411	IF(n0.gt.0) THEN
[524]	412	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
	413	ccc & ,'contenu de la maille : ',n0
	414
	415	DO l=1,llm
	416	IF(nl(l).gt.0) THEN
	417	iju=0
	418	c indicage des mailles concernees par le traitement special
	419	DO ij=iip2,ip1jm
	420	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	421	iju=iju+1
	422	indu(iju)=ij
	423	ENDIF
	424	ENDDO
	425	niju=iju
	426	c PRINT*,'niju,nl',niju,nl(l)
	427
	428	c traitement des mailles
	429	DO iju=1,niju
	430	ij=indu(iju)
	431	j=(ij-1)/iip1+1
	432	zu_m=u_m(ij,l)
	433	u_mq(ij,l)=0.
	434	IF(zu_m.gt.0.) THEN
	435	ijq=ij
	436	i=ijq-(j-1)*iip1
	437	c accumulation pour les mailles completements advectees
[2270]	438	do while(zu_m.gt.masse(ijq,l,iq))
	439	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
	440	& *masse(ijq,l,iq)
	441	zu_m=zu_m-masse(ijq,l,iq)
[524]	442	i=mod(i-2+iim,iim)+1
	443	ijq=(j-1)*iip1+i
	444	ENDDO
	445	c ajout de la maille non completement advectee
	446	u_mq(ij,l)=u_mq(ij,l)+zu_m*
[2270]	447	& (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
	448	& *dxq(ijq,l))
[524]	449	ELSE
	450	ijq=ij+1
	451	i=ijq-(j-1)*iip1
	452	c accumulation pour les mailles completements advectees
[2270]	453	do while(-zu_m.gt.masse(ijq,l,iq))
	454	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
	455	& *masse(ijq,l,iq)
	456	zu_m=zu_m+masse(ijq,l,iq)
[524]	457	i=mod(i,iim)+1
	458	ijq=(j-1)*iip1+i
	459	ENDDO
	460	c ajout de la maille non completement advectee
[2270]	461	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
	462	& 0.5(1.+zu_m/masse(ijq,l,iq))dxq(ijq,l))
[524]	463	ENDIF
	464	ENDDO
	465	ENDIF
	466	ENDDO
	467	ENDIF ! n0.gt.0
	468
	469
	470
	471	c bouclage en latitude
	472
	473	DO l=1,llm
	474	DO ij=iip1+iip1,ip1jm,iip1
	475	u_mq(ij,l)=u_mq(ij-iim,l)
	476	ENDDO
	477	ENDDO
	478
[2270]	479	! CRisi: appel récursif de l'advection sur les fils.
	480	! Il faut faire ça avant d'avoir mis à jour q et masse
[4368]	481	!write(,) 'vlspltqs 326: iq,nqChildren(iq)=',iq,
	482	! & tracers(iq)%nqChildren
[2270]	483
[4368]	484	do ifils=1,tracers(iq)%nqDescen
	485	iq2=tracers(iq)%iqDescen(ifils)
	486	DO l=1,llm
[2270]	487	DO ij=iip2,ip1jm
[4368]	488	! On a besoin de q et masse seulement entre iip2 et ip1jm
	489	masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
	490	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
[2270]	491	enddo
[4368]	492	enddo
	493	enddo
	494	do ifils=1,tracers(iq)%nqChildren
	495	iq2=tracers(iq)%iqDescen(ifils)
	496	call vlx(Ratio,pente_max,masseq,u_mq,iq2)
	497	enddo
[2270]	498	! end CRisi
[524]	499
	500	c calcul des tendances
	501
	502	DO l=1,llm
	503	DO ij=iip2+1,ip1jm
[2270]	504	new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
	505	q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
[524]	506	& u_mq(ij-1,l)-u_mq(ij,l))
	507	& /new_m
[2270]	508	masse(ij,l,iq)=new_m
[524]	509	ENDDO
	510	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	511	DO ij=iip1+iip1,ip1jm,iip1
[2270]	512	q(ij-iim,l,iq)=q(ij,l,iq)
	513	masse(ij-iim,l,iq)=masse(ij,l,iq)
[524]	514	ENDDO
	515	ENDDO
	516
[2270]	517	! retablir les fils en rapport de melange par rapport a l'air:
	518	! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
	519	! puis on boucle en longitude
[4368]	520	do ifils=1,tracers(iq)%nqDescen
	521	iq2=tracers(iq)%iqDescen(ifils)
	522	DO l=1,llm
[2270]	523	DO ij=iip2+1,ip1jm
	524	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
	525	enddo
	526	DO ij=iip1+iip1,ip1jm,iip1
[4368]	527	q(ij-iim,l,iq2)=q(ij,l,iq2)
	528	enddo
	529	enddo
	530	enddo
[2270]	531
[524]	532	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	533	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	534
	535
	536	RETURN
	537	END
[2270]	538	SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
[4368]	539	USE infotrac, ONLY : nqtot,tracers ! CRisi
[524]	540	c
	541	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	542	c
	543	c ********************************************************************
	544	c Shema d'advection " pseudo amont " .
	545	c ********************************************************************
	546	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	547	c qsat est un argument de sortie pour le s-pg ....
	548	c
	549	c
	550	c --------------------------------------------------------------------
[2597]	551
	552	USE comconst_mod, ONLY: pi
	553
[524]	554	IMPLICIT NONE
	555	c
[2597]	556	include "dimensions.h"
	557	include "paramet.h"
	558	include "comgeom.h"
[524]	559	c
	560	c
	561	c Arguments:
	562	c ----------
[2270]	563	REAL masse(ip1jmp1,llm,nqtot),pente_max
[524]	564	REAL masse_adv_v( ip1jm,llm)
[2270]	565	REAL q(ip1jmp1,llm,nqtot)
[524]	566	REAL qsat(ip1jmp1,llm)
[2270]	567	INTEGER iq ! CRisi
[524]	568	c
	569	c Local
	570	c ---------
	571	c
	572	INTEGER i,ij,l
	573	c
	574	REAL airej2,airejjm,airescb(iim),airesch(iim)
	575	REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
	576	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
	577	REAL qbyv(ip1jm,llm)
	578
	579	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	580	c REAL newq,oldmasse
	581	Logical first,testcpu
	582	REAL temps0,temps1,temps2,temps3,temps4,temps5
	583	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	584	SAVE first,testcpu
	585
	586	REAL convpn,convps,convmpn,convmps
	587	REAL sinlon(iip1),sinlondlon(iip1)
	588	REAL coslon(iip1),coslondlon(iip1)
	589	SAVE sinlon,coslon,sinlondlon,coslondlon
	590	SAVE airej2,airejjm
[2270]	591
	592	REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
	593	INTEGER ifils,iq2 ! CRisi
[524]	594	c
	595	c
	596	REAL SSUM
	597
	598	DATA first,testcpu/.true.,.false./
	599	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
	600
	601	IF(first) THEN
	602	PRINT*,'Shema Amont nouveau appele dans Vanleer '
	603	first=.false.
	604	do i=2,iip1
	605	coslon(i)=cos(rlonv(i))
	606	sinlon(i)=sin(rlonv(i))
	607	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	608	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	609	ENDDO
	610	coslon(1)=coslon(iip1)
	611	coslondlon(1)=coslondlon(iip1)
	612	sinlon(1)=sinlon(iip1)
	613	sinlondlon(1)=sinlondlon(iip1)
	614	airej2 = SSUM( iim, aire(iip2), 1 )
	615	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
	616	ENDIF
	617
	618	c
	619
	620
	621	DO l = 1, llm
	622	c
	623	c --------------------------------
	624	c CALCUL EN LATITUDE
	625	c --------------------------------
	626
	627	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
	628	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
	629	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
	630
	631	DO i = 1, iim
[2270]	632	airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
	633	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
[524]	634	ENDDO
	635	qpns = SSUM( iim, airescb ,1 ) / airej2
	636	qpsn = SSUM( iim, airesch ,1 ) / airejjm
	637
	638	c calcul des pentes aux points v
	639
	640	DO ij=1,ip1jm
[2270]	641	dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
[524]	642	adyqv(ij)=abs(dyqv(ij))
	643	ENDDO
	644
	645	c calcul des pentes aux points scalaires
	646
	647	DO ij=iip2,ip1jm
	648	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
	649	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
	650	dyqmax(ij)=pente_max*dyqmax(ij)
	651	ENDDO
	652
	653	c calcul des pentes aux poles
	654
	655	DO ij=1,iip1
[2270]	656	dyq(ij,l)=qpns-q(ij+iip1,l,iq)
	657	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
[524]	658	ENDDO
	659
	660	c filtrage de la derivee
	661	dyn1=0.
	662	dys1=0.
	663	dyn2=0.
	664	dys2=0.
	665	DO ij=1,iim
	666	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
	667	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
	668	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
	669	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
	670	ENDDO
	671	DO ij=1,iip1
	672	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
	673	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
	674	ENDDO
	675
	676	c calcul des pentes limites aux poles
	677
	678	fn=1.
	679	fs=1.
	680	DO ij=1,iim
	681	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
	682	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
	683	ENDIF
	684	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
	685	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
	686	ENDIF
	687	ENDDO
	688	DO ij=1,iip1
	689	dyq(ij,l)=fn*dyq(ij,l)
	690	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
	691	ENDDO
	692
	693	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	694	C En memoire de dIFferents tests sur la
	695	C limitation des pentes aux poles.
	696	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	697	C PRINT*,dyq(1)
	698	C PRINT*,dyqv(iip1+1)
[1520]	699	C appn=abs(dyq(1)/dyqv(iip1+1))
[524]	700	C PRINT*,dyq(ip1jm+1)
	701	C PRINT*,dyqv(ip1jm-iip1+1)
[1520]	702	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[524]	703	C DO ij=2,iim
[1520]	704	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
	705	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
[524]	706	C ENDDO
[1520]	707	C appn=min(pente_max/appn,1.)
	708	C apps=min(pente_max/apps,1.)
[524]	709	C
	710	C
	711	C cas ou on a un extremum au pole
	712	C
	713	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
[1520]	714	C & appn=0.
[524]	715	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	716	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
[1520]	717	C & apps=0.
[524]	718	C
	719	C limitation des pentes aux poles
	720	C DO ij=1,iip1
[1520]	721	C dyq(ij)=appn*dyq(ij)
	722	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
[524]	723	C ENDDO
	724	C
	725	C test
	726	C DO ij=1,iip1
	727	C dyq(iip1+ij)=0.
	728	C dyq(ip1jm+ij-iip1)=0.
	729	C ENDDO
	730	C DO ij=1,ip1jmp1
	731	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
	732	C ENDDO
	733	C
	734	C changement 10 07 96
	735	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	736	C & THEN
	737	C DO ij=1,iip1
	738	C dyqmax(ij)=0.
	739	C ENDDO
	740	C ELSE
	741	C DO ij=1,iip1
	742	C dyqmax(ij)=pente_max*abs(dyqv(ij))
	743	C ENDDO
	744	C ENDIF
	745	C
	746	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	747	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	748	C &THEN
	749	C DO ij=ip1jm+1,ip1jmp1
	750	C dyqmax(ij)=0.
	751	C ENDDO
	752	C ELSE
	753	C DO ij=ip1jm+1,ip1jmp1
	754	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
	755	C ENDDO
	756	C ENDIF
	757	C fin changement 10 07 96
	758	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	759
	760	c calcul des pentes limitees
	761
	762	DO ij=iip2,ip1jm
	763	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
	764	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
	765	ELSE
	766	dyq(ij,l)=0.
	767	ENDIF
	768	ENDDO
	769
	770	ENDDO
	771
	772	DO l=1,llm
	773	DO ij=1,ip1jm
	774	IF( masse_adv_v(ij,l).GT.0. ) THEN
[2270]	775	qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq ) +
	776	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)
	777	, /masse(ij+iip1,l,iq)))
[524]	778	ELSE
[2270]	779	qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
	780	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
[524]	781	ENDIF
	782	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
	783	ENDDO
	784	ENDDO
	785
	786
[2270]	787	! CRisi: appel récursif de l'advection sur les fils.
	788	! Il faut faire ça avant d'avoir mis à jour q et masse
[4368]	789	!write(,) 'vlyqs 689: iq,nqChildren(iq)=',iq,
	790	! & tracers(iq)%nqChildren
[2270]	791
[4368]	792	do ifils=1,tracers(iq)%nqDescen
	793	iq2=tracers(iq)%iqDescen(ifils)
	794	DO l=1,llm
	795	DO ij=1,ip1jmp1
	796	masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
	797	Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
[2270]	798	enddo
[4368]	799	enddo
	800	enddo
	801	do ifils=1,tracers(iq)%nqChildren
	802	iq2=tracers(iq)%iqDescen(ifils)
	803	!write(,) 'vlyqs 783: appel rec de vly, iq2=',iq2
	804	call vly(Ratio,pente_max,masseq,qbyv,iq2)
	805	enddo
[2270]	806
[524]	807	DO l=1,llm
	808	DO ij=iip2,ip1jm
[2270]	809	newmasse=masse(ij,l,iq)
[524]	810	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
[2270]	811	q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
	812	& -qbyv(ij-iip1,l))/newmasse
	813	masse(ij,l,iq)=newmasse
[524]	814	ENDDO
	815	c.-. ancienne version
	816	convpn=SSUM(iim,qbyv(1,l),1)/apoln
	817	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
	818	DO ij = 1,iip1
[2270]	819	newmasse=masse(ij,l,iq)+convmpn*aire(ij)
	820	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpnaire(ij))/
[524]	821	& newmasse
[2270]	822	masse(ij,l,iq)=newmasse
[524]	823	ENDDO
	824	convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
	825	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
	826	DO ij = ip1jm+1,ip1jmp1
[2270]	827	newmasse=masse(ij,l,iq)+convmps*aire(ij)
	828	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpsaire(ij))/
[524]	829	& newmasse
[2270]	830	masse(ij,l,iq)=newmasse
[524]	831	ENDDO
	832	c.-. fin ancienne version
	833
	834	c._. nouvelle version
	835	c convpn=SSUM(iim,qbyv(1,l),1)
	836	c convmpn=ssum(iim,masse_adv_v(1,l),1)
	837	c oldmasse=ssum(iim,masse(1,l),1)
	838	c newmasse=oldmasse+convmpn
	839	c newq=(q(1,l)*oldmasse+convpn)/newmasse
	840	c newmasse=newmasse/apoln
	841	c DO ij = 1,iip1
	842	c q(ij,l)=newq
[2270]	843	c masse(ij,l,iq)=newmasse*aire(ij)
[524]	844	c ENDDO
	845	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	846	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	847	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
	848	c newmasse=oldmasse+convmps
	849	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
	850	c newmasse=newmasse/apols
	851	c DO ij = ip1jm+1,ip1jmp1
	852	c q(ij,l)=newq
[2270]	853	c masse(ij,l,iq)=newmasse*aire(ij)
[524]	854	c ENDDO
	855	c._. fin nouvelle version
	856	ENDDO
	857
[2286]	858	!write(,) 'vly 866'
[2270]	859
	860	! retablir les fils en rapport de melange par rapport a l'air:
[4368]	861	do ifils=1,tracers(iq)%nqDescen
	862	iq2=tracers(iq)%iqDescen(ifils)
	863	DO l=1,llm
[2270]	864	DO ij=1,ip1jmp1
	865	q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
	866	enddo
[4368]	867	enddo
	868	enddo
[2286]	869	!write(,) 'vly 879'
[2270]	870
[524]	871	RETURN
	872	END

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

Download in other formats:

Original Format