Context Navigation

vlspltqs.F @ 5394

Last change on this file since 5394 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
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: 25.0 KB

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

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

Download in other formats:

Original Format