Context Navigation

vlspltqs.F @ 5441

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

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

Download in other formats:

Original Format