Context Navigation

vlspltqs_p.F @ 4289

Last change on this file since 4289 was 709, checked in by Laurent Fairhead, 18 years ago
Nouvelles versions de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 24.6 KB

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

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

Download in other formats:

Original Format