Context Navigation

vlspltqs_loc.F @ 5456

Last change on this file since 5456 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
File size: 23.3 KB

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

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

Download in other formats:

Original Format