Context Navigation

vlsplt_p.F @ 839

Last change on this file since 839 was 764, checked in by Laurent Fairhead, 17 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.8 KB

Rev	Line
[630]	1	SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
	2	c
	3	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	4	c
	5	c ********************************************************************
	6	c Shema d'advection " pseudo amont " .
	7	c ********************************************************************
	8	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	9	c
	10	c pente_max facteur de limitation des pentes: 2 en general
	11	c 0 pour un schema amont
	12	c pbaru,pbarv,w flux de masse en u ,v ,w
	13	c pdt pas de temps
	14	c
	15	c --------------------------------------------------------------------
	16	USE parallel
	17	USE mod_hallo
	18	USE Vampir
	19	IMPLICIT NONE
	20	c
	21	#include "dimensions.h"
	22	#include "paramet.h"
	23	#include "logic.h"
	24	#include "comvert.h"
	25	#include "comconst.h"
	26
	27	c
	28	c Arguments:
	29	c ----------
	30	REAL masse(ip1jmp1,llm),pente_max
	31	c REAL masse(iip1,jjp1,llm),pente_max
	32	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
	33	REAL q(ip1jmp1,llm)
	34	c REAL q(iip1,jjp1,llm)
	35	REAL w(ip1jmp1,llm),pdt
	36	c
	37	c Local
	38	c ---------
	39	c
	40	INTEGER i,ij,l,j,ii
	41	INTEGER ijlqmin,iqmin,jqmin,lqmin
	42	c
	43	REAL zm(ip1jmp1,llm),newmasse
	44	REAL mu(ip1jmp1,llm)
	45	REAL mv(ip1jm,llm)
	46	REAL mw(ip1jmp1,llm+1)
	47	REAL zq(ip1jmp1,llm),zz
	48	REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
	49	REAL second,temps0,temps1,temps2,temps3
	50	REAL ztemps1,ztemps2,ztemps3
	51	REAL zzpbar, zzw
	52	LOGICAL testcpu
	53	SAVE testcpu
	54	SAVE temps1,temps2,temps3
	55	INTEGER iminn,imaxx
	56
	57	REAL qmin,qmax
	58	DATA qmin,qmax/0.,1.e33/
	59	DATA testcpu/.false./
	60	DATA temps1,temps2,temps3/0.,0.,0./
	61	INTEGER ijb,ije
	62	type(request) :: MyRequest1
	63	type(request) :: MyRequest2
	64
	65	call SetTag(MyRequest1,100)
	66	call SetTag(MyRequest2,101)
	67
	68	zzpbar = 0.5 * pdt
	69	zzw = pdt
	70
	71	ijb=ij_begin
	72	ije=ij_end
	73	if (pole_nord) ijb=ijb+iip1
	74	if (pole_sud) ije=ije-iip1
	75
	76	DO l=1,llm
	77	DO ij = ijb,ije
	78	mu(ij,l)=pbaru(ij,l) * zzpbar
	79	ENDDO
	80	ENDDO
	81
	82	ijb=ij_begin-iip1
	83	ije=ij_end
	84	if (pole_nord) ijb=ij_begin
	85	if (pole_sud) ije=ij_end-iip1
	86
	87	DO l=1,llm
	88	DO ij=ijb,ije
	89	mv(ij,l)=pbarv(ij,l) * zzpbar
	90	ENDDO
	91	ENDDO
	92
	93	ijb=ij_begin
	94	ije=ij_end
	95
	96	DO l=1,llm
	97	DO ij=ijb,ije
	98	mw(ij,l)=w(ij,l) * zzw
	99	ENDDO
	100	ENDDO
	101
	102	DO ij=ijb,ije
	103	mw(ij,llm+1)=0.
	104	ENDDO
	105
	106	c CALL SCOPY(ijp1llm,q,1,zq,1)
	107	c CALL SCOPY(ijp1llm,masse,1,zm,1)
	108
	109	ijb=ij_begin
	110	ije=ij_end
	111	zq(ijb:ije,:)=q(ijb:ije,:)
	112	zm(ijb:ije,:)=masse(ijb:ije,:)
	113
	114
	115	c print*,'Entree vlx1'
	116	c call minmaxq(zq,qmin,qmax,'avant vlx ')
	117	call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1)
	118	call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end)
	119	call VTb(VTHallo)
	120	call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
	121	call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest1)
	122	call SendRequest(MyRequest1)
	123	call VTe(VTHallo)
	124	call vlx_p(zq,pente_max,zm,mu,ij_begin+2iip1,ij_end-2iip1)
	125	c call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
	126	call VTb(VTHallo)
	127	call WaitRecvRequest(MyRequest1)
	128	call VTe(VTHallo)
	129
	130
	131	c print*,'Sortie vlx1'
	132	c call minmaxq(zq,qmin,qmax,'apres vlx1 ')
	133
	134	c print*,'Entree vly1'
	135	c call exchange_hallo(zq,ip1jmp1,llm,2,2)
	136	c call exchange_hallo(zm,ip1jmp1,llm,1,1)
	137
	138	call vly_p(zq,pente_max,zm,mv)
	139	c call minmaxq(zq,qmin,qmax,'apres vly1 ')
	140	c print*,'Sortie vly1'
	141	call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
	142	call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
	143	call VTb(VTHallo)
	144	call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
	145	call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest2)
	146	call SendRequest(MyRequest2)
	147	call VTe(VTHallo)
	148	call vlz_p(zq,pente_max,zm,mw,ij_begin+2iip1,ij_end-2iip1)
	149	call VTb(VTHallo)
	150	call WaitRecvRequest(MyRequest2)
	151
	152	call VTe(VTHallo)
	153
	154	c call minmaxq(zq,qmin,qmax,'apres vlz ')
	155
	156
	157
	158
	159	call vly_p(zq,pente_max,zm,mv)
	160	c call minmaxq(zq,qmin,qmax,'apres vly ')
	161
	162
	163	call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
	164	c call minmaxq(zq,qmin,qmax,'apres vlx2 ')
	165
	166
	167	ijb=ij_begin
	168	ije=ij_end
	169
	170	DO l=1,llm
	171	DO ij=ijb,ije
	172	q(ij,l)=zq(ij,l)
	173	ENDDO
	174	ENDDO
	175
	176
	177	DO l=1,llm
	178	DO ij=ijb,ije-iip1+1,iip1
	179	q(ij+iim,l)=q(ij,l)
	180	ENDDO
	181	ENDDO
	182
	183	call WaitSendRequest(MyRequest1)
	184	call WaitSendRequest(MyRequest2)
	185
	186	RETURN
	187	END
	188
	189
	190	SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x)
	191
	192	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	193	c
	194	c ********************************************************************
	195	c Shema d'advection " pseudo amont " .
	196	c ********************************************************************
	197	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	198	c
	199	c
	200	c --------------------------------------------------------------------
	201	USE Parallel
	202	IMPLICIT NONE
	203	c
	204	#include "dimensions.h"
	205	#include "paramet.h"
	206	#include "logic.h"
	207	#include "comvert.h"
	208	#include "comconst.h"
	209	c
	210	c
	211	c Arguments:
	212	c ----------
	213	REAL masse(ip1jmp1,llm),pente_max
	214	REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
	215	REAL q(ip1jmp1,llm)
	216	REAL w(ip1jmp1,llm)
	217	c
	218	c Local
	219	c ---------
	220	c
	221	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
	222	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
	223	c
	224	REAL new_m,zu_m,zdum(ip1jmp1,llm)
	225	REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
	226	REAL zz(ip1jmp1)
	227	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
	228	REAL u_mq(ip1jmp1,llm)
	229
[764]	230	Logical extremum
[630]	231
	232	REAL SSUM
	233	EXTERNAL SSUM
	234
	235	REAL z1,z2,z3
	236
	237	INTEGER ijb,ije,ijb_x,ije_x
	238
	239	c calcul de la pente a droite et a gauche de la maille
	240
	241	ijb=ijb_x
	242	ije=ije_x
	243
	244	if (pole_nord.and.ijb==1) ijb=ijb+iip1
	245	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
	246
	247	IF (pente_max.gt.-1.e-5) THEN
	248	c IF (pente_max.gt.10) THEN
	249
	250	c calcul des pentes avec limitation, Van Leer scheme I:
	251	c -----------------------------------------------------
	252
	253	c calcul de la pente aux points u
[764]	254	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	255	DO l = 1, llm
	256
	257	DO ij=ijb,ije-1
	258	dxqu(ij)=q(ij+1,l)-q(ij,l)
	259	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	260	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	261	ENDDO
	262	DO ij=ijb+iip1-1,ije,iip1
	263	dxqu(ij)=dxqu(ij-iim)
	264	c sigu(ij)=sigu(ij-iim)
	265	ENDDO
	266
	267	DO ij=ijb,ije
	268	adxqu(ij)=abs(dxqu(ij))
	269	ENDDO
	270
	271	c calcul de la pente maximum dans la maille en valeur absolue
	272
	273	DO ij=ijb+1,ije
	274	dxqmax(ij,l)=pente_max*
	275	, min(adxqu(ij-1),adxqu(ij))
	276	c limitation subtile
	277	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	278
	279
	280	ENDDO
	281
	282	DO ij=ijb+iip1-1,ije,iip1
	283	dxqmax(ij-iim,l)=dxqmax(ij,l)
	284	ENDDO
	285
	286	DO ij=ijb+1,ije
	287	#ifdef CRAY
	288	dxq(ij,l)=
	289	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	290	#else
	291	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	292	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	293	ELSE
	294	c extremum local
	295	dxq(ij,l)=0.
	296	ENDIF
	297	#endif
	298	dxq(ij,l)=0.5*dxq(ij,l)
	299	dxq(ij,l)=
	300	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	301	ENDDO
	302
	303	ENDDO ! l=1,llm
[764]	304	c$OMP END DO NOWAIT
[630]	305	c print*,'Ok calcul des pentes'
	306
	307	ELSE ! (pente_max.lt.-1.e-5)
	308
	309	c Pentes produits:
	310	c ----------------
[764]	311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	312	DO l = 1, llm
	313	DO ij=ijb,ije-1
	314	dxqu(ij)=q(ij+1,l)-q(ij,l)
	315	ENDDO
	316	DO ij=ijb+iip1-1,ije,iip1
	317	dxqu(ij)=dxqu(ij-iim)
	318	ENDDO
	319
	320	DO ij=ijb+1,ije
	321	zz(ij)=dxqu(ij-1)*dxqu(ij)
	322	zz(ij)=zz(ij)+zz(ij)
	323	IF(zz(ij).gt.0) THEN
	324	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	325	ELSE
	326	c extremum local
	327	dxq(ij,l)=0.
	328	ENDIF
	329	ENDDO
	330
	331	ENDDO
[764]	332	c$OMP END DO NOWAIT
[630]	333	ENDIF ! (pente_max.lt.-1.e-5)
	334
	335	c bouclage de la pente en iip1:
	336	c -----------------------------
[764]	337	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	338	DO l=1,llm
	339	DO ij=ijb+iip1-1,ije,iip1
	340	dxq(ij-iim,l)=dxq(ij,l)
	341	ENDDO
[764]	342	DO ij=ijb,ije
[630]	343	iadvplus(ij,l)=0
	344	ENDDO
	345
	346	ENDDO
[764]	347	c$OMP END DO NOWAIT
[630]	348	c print*,'Bouclage en iip1'
	349
	350	c calcul des flux a gauche et a droite
	351
	352	#ifdef CRAY
[764]	353	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	354	DO l=1,llm
	355	DO ij=ijb,ije-1
	356	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	357	, 1.+u_m(ij,l)/masse(ij+1,l),
	358	, u_m(ij,l))
	359	zdum(ij,l)=0.5*zdum(ij,l)
	360	u_mq(ij,l)=cvmgp(
	361	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	362	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	363	, u_m(ij,l))
	364	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	365	ENDDO
	366	ENDDO
[764]	367	c$OMP END DO NOWAIT
[630]	368	#else
	369	c on cumule le flux correspondant a toutes les mailles dont la masse
	370	c au travers de la paroi pENDant le pas de temps.
	371	c print*,'Cumule ....'
[764]	372	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	373	DO l=1,llm
	374	DO ij=ijb,ije-1
	375	c print*,'masse(',ij,')=',masse(ij,l)
	376	IF (u_m(ij,l).gt.0.) THEN
	377	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	378	u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
	379	ELSE
	380	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	381	u_mq(ij,l)=u_m(ij,l)(q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
	382	ENDIF
	383	ENDDO
	384	ENDDO
[764]	385	c$OMP END DO NOWAIT
[630]	386	#endif
	387	c stop
	388
	389	c go to 9999
	390	c detection des points ou on advecte plus que la masse de la
	391	c maille
[764]	392	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	393	DO l=1,llm
	394	DO ij=ijb,ije-1
	395	IF(zdum(ij,l).lt.0) THEN
	396	iadvplus(ij,l)=1
	397	u_mq(ij,l)=0.
	398	ENDIF
	399	ENDDO
	400	ENDDO
[764]	401	c$OMP END DO NOWAIT
[630]	402	c print*,'Ok test 1'
[764]	403	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	404	DO l=1,llm
	405	DO ij=ijb+iip1-1,ije,iip1
	406	iadvplus(ij,l)=iadvplus(ij-iim,l)
	407	ENDDO
	408	ENDDO
[764]	409	c$OMP END DO NOWAIT
[630]	410	c print*,'Ok test 2'
	411
	412
	413	c traitement special pour le cas ou on advecte en longitude plus que le
	414	c contenu de la maille.
	415	c cette partie est mal vectorisee.
	416
	417	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	418
	419	n0=0
[764]	420	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	421	DO l=1,llm
	422	nl(l)=0
	423	DO ij=ijb,ije
	424	nl(l)=nl(l)+iadvplus(ij,l)
	425	ENDDO
	426	n0=n0+nl(l)
	427	ENDDO
[764]	428	c$OMP END DO NOWAIT
[630]	429	cym IF(n0.gt.1) THEN
[764]	430	cym IF(n0.gt.0) THEN
[630]	431
	432	c PRINT*,'Nombre de points pour lesquels on advect plus que le'
	433	c & ,'contenu de la maille : ',n0
[764]	434	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	435	DO l=1,llm
	436	IF(nl(l).gt.0) THEN
	437	iju=0
	438	c indicage des mailles concernees par le traitement special
	439	DO ij=ijb,ije
	440	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	441	iju=iju+1
	442	indu(iju)=ij
	443	ENDIF
	444	ENDDO
	445	niju=iju
	446	c PRINT*,'niju,nl',niju,nl(l)
	447
	448	c traitement des mailles
	449	DO iju=1,niju
	450	ij=indu(iju)
	451	j=(ij-1)/iip1+1
	452	zu_m=u_m(ij,l)
	453	u_mq(ij,l)=0.
	454	IF(zu_m.gt.0.) THEN
	455	ijq=ij
	456	i=ijq-(j-1)*iip1
	457	c accumulation pour les mailles completements advectees
	458	do while(zu_m.gt.masse(ijq,l))
	459	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	460	zu_m=zu_m-masse(ijq,l)
	461	i=mod(i-2+iim,iim)+1
	462	ijq=(j-1)*iip1+i
	463	ENDDO
	464	c ajout de la maille non completement advectee
	465	u_mq(ij,l)=u_mq(ij,l)+zu_m*
	466	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
	467	ELSE
	468	ijq=ij+1
	469	i=ijq-(j-1)*iip1
	470	c accumulation pour les mailles completements advectees
	471	do while(-zu_m.gt.masse(ijq,l))
	472	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	473	zu_m=zu_m+masse(ijq,l)
	474	i=mod(i,iim)+1
	475	ijq=(j-1)*iip1+i
	476	ENDDO
	477	c ajout de la maille non completement advectee
	478	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
	479	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
	480	ENDIF
	481	ENDDO
	482	ENDIF
	483	ENDDO
[764]	484	c$OMP END DO NOWAIT
	485	cym ENDIF ! n0.gt.0
[630]	486	9999 continue
	487
	488
	489	c bouclage en latitude
	490	c print*,'Avant bouclage en latitude'
[764]	491	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	492	DO l=1,llm
	493	DO ij=ijb+iip1-1,ije,iip1
	494	u_mq(ij,l)=u_mq(ij-iim,l)
	495	ENDDO
	496	ENDDO
[764]	497	c$OMP END DO NOWAIT
[630]	498
	499	c calcul des tENDances
[764]	500	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	501	DO l=1,llm
	502	DO ij=ijb+1,ije
	503	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	504	q(ij,l)=(q(ij,l)*masse(ij,l)+
	505	& u_mq(ij-1,l)-u_mq(ij,l))
	506	& /new_m
	507	masse(ij,l)=new_m
	508	ENDDO
	509	c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	510	DO ij=ijb+iip1-1,ije,iip1
	511	q(ij-iim,l)=q(ij,l)
	512	masse(ij-iim,l)=masse(ij,l)
	513	ENDDO
	514	ENDDO
[764]	515	c$OMP END DO NOWAIT
[630]	516	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	517	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	518
	519
	520	RETURN
	521	END
	522
	523
	524	SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
	525	c
	526	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	527	c
	528	c ********************************************************************
	529	c Shema d'advection " pseudo amont " .
	530	c ********************************************************************
	531	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	532	c dq sont des arguments de sortie pour le s-pg ....
	533	c
	534	c
	535	c --------------------------------------------------------------------
	536	USE parallel
	537	IMPLICIT NONE
	538	c
	539	#include "dimensions.h"
	540	#include "paramet.h"
	541	#include "logic.h"
	542	#include "comvert.h"
	543	#include "comconst.h"
	544	#include "comgeom.h"
	545	c
	546	c
	547	c Arguments:
	548	c ----------
	549	REAL masse(ip1jmp1,llm),pente_max
	550	REAL masse_adv_v( ip1jm,llm)
	551	REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
	552	c
	553	c Local
	554	c ---------
	555	c
	556	INTEGER i,ij,l
	557	c
	558	REAL airej2,airejjm,airescb(iim),airesch(iim)
	559	REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
	560	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
	561	REAL qbyv(ip1jm,llm)
	562
	563	REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	564	c REAL newq,oldmasse
	565	Logical extremum,first,testcpu
	566	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
	567	SAVE temps0,temps1,temps2,temps3,temps4,temps5
[764]	568	c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
[630]	569	SAVE first,testcpu
[764]	570	c$OMP THREADPRIVATE(first,testcpu)
[630]	571
	572	REAL convpn,convps,convmpn,convmps
	573	real massepn,masseps,qpn,qps
	574	REAL sinlon(iip1),sinlondlon(iip1)
	575	REAL coslon(iip1),coslondlon(iip1)
	576	SAVE sinlon,coslon,sinlondlon,coslondlon
[764]	577	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
[630]	578	SAVE airej2,airejjm
[764]	579	c$OMP THREADPRIVATE(airej2,airejjm)
[630]	580	c
	581	c
	582	REAL SSUM
	583	EXTERNAL SSUM
	584
	585	DATA first,testcpu/.true.,.false./
	586	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
	587	INTEGER ijb,ije
	588
	589	IF(first) THEN
	590	c PRINT*,'Shema Amont nouveau appele dans Vanleer '
	591	first=.false.
	592	do i=2,iip1
	593	coslon(i)=cos(rlonv(i))
	594	sinlon(i)=sin(rlonv(i))
	595	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	596	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	597	ENDDO
	598	coslon(1)=coslon(iip1)
	599	coslondlon(1)=coslondlon(iip1)
	600	sinlon(1)=sinlon(iip1)
	601	sinlondlon(1)=sinlondlon(iip1)
	602	airej2 = SSUM( iim, aire(iip2), 1 )
	603	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
	604	ENDIF
	605
	606	c
	607	c PRINT*,'CALCUL EN LATITUDE'
	608
[764]	609	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	610	DO l = 1, llm
	611	c
	612	c --------------------------------
	613	c CALCUL EN LATITUDE
	614	c --------------------------------
	615
	616	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
	617	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
	618	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
	619
	620	if (pole_nord) then
	621	DO i = 1, iim
	622	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
	623	ENDDO
	624	qpns = SSUM( iim, airescb ,1 ) / airej2
	625	endif
	626
	627	if (pole_sud) then
	628	DO i = 1, iim
	629	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
	630	ENDDO
	631	qpsn = SSUM( iim, airesch ,1 ) / airejjm
	632	endif
	633
	634
	635
	636	c calcul des pentes aux points v
	637
	638	ijb=ij_begin-2*iip1
	639	ije=ij_end+iip1
	640	if (pole_nord) ijb=ij_begin
	641	if (pole_sud) ije=ij_end-iip1
	642
	643	DO ij=ijb,ije
	644	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
	645	adyqv(ij)=abs(dyqv(ij))
	646	ENDDO
	647
	648	c calcul des pentes aux points scalaires
	649	ijb=ij_begin-iip1
	650	ije=ij_end+iip1
	651	if (pole_nord) ijb=ij_begin+iip1
	652	if (pole_sud) ije=ij_end-iip1
	653
	654	DO ij=ijb,ije
	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	IF (pole_nord) THEN
	662	DO ij=1,iip1
	663	dyq(ij,l)=qpns-q(ij+iip1,l)
	664	ENDDO
	665
	666	dyn1=0.
	667	dyn2=0.
	668	DO ij=1,iim
	669	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
	670	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
	671	ENDDO
	672	DO ij=1,iip1
	673	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
	674	ENDDO
	675
	676	DO ij=1,iip1
	677	dyq(ij,l)=0.
	678	ENDDO
	679	c ym tout cela ne sert pas a grand chose
	680	ENDIF
	681
	682	IF (pole_sud) THEN
	683
	684	DO ij=1,iip1
	685	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
	686	ENDDO
	687
	688	dys1=0.
	689	dys2=0.
	690
	691	DO ij=1,iim
	692	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
	693	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
	694	ENDDO
	695
	696	DO ij=1,iip1
	697	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
	698	ENDDO
	699
	700	DO ij=1,iip1
	701	dyq(ip1jm+ij,l)=0.
	702	ENDDO
	703	c ym tout cela ne sert pas a grand chose
	704	ENDIF
	705
	706	c filtrage de la derivee
	707
	708	c calcul des pentes limites aux poles
	709	c ym partie inutile
	710	c goto 8888
	711	c fn=1.
	712	c fs=1.
	713	c DO ij=1,iim
	714	c IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
	715	c fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
	716	c ENDIF
	717	c IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
	718	c fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
	719	c ENDIF
	720	c ENDDO
	721	c DO ij=1,iip1
	722	c dyq(ij,l)=fn*dyq(ij,l)
	723	c dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
	724	c ENDDO
	725	c 8888 continue
	726
	727
	728	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	729	C En memoire de dIFferents tests sur la
	730	C limitation des pentes aux poles.
	731	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	732	C PRINT*,dyq(1)
	733	C PRINT*,dyqv(iip1+1)
	734	C apn=abs(dyq(1)/dyqv(iip1+1))
	735	C PRINT*,dyq(ip1jm+1)
	736	C PRINT*,dyqv(ip1jm-iip1+1)
	737	C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
	738	C DO ij=2,iim
	739	C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
	740	C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
	741	C ENDDO
	742	C apn=min(pente_max/apn,1.)
	743	C aps=min(pente_max/aps,1.)
	744	C
	745	C
	746	C cas ou on a un extremum au pole
	747	C
	748	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	749	C & apn=0.
	750	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	751	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	752	C & aps=0.
	753	C
	754	C limitation des pentes aux poles
	755	C DO ij=1,iip1
	756	C dyq(ij)=apn*dyq(ij)
	757	C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
	758	C ENDDO
	759	C
	760	C test
	761	C DO ij=1,iip1
	762	C dyq(iip1+ij)=0.
	763	C dyq(ip1jm+ij-iip1)=0.
	764	C ENDDO
	765	C DO ij=1,ip1jmp1
	766	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
	767	C ENDDO
	768	C
	769	C changement 10 07 96
	770	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	771	C & THEN
	772	C DO ij=1,iip1
	773	C dyqmax(ij)=0.
	774	C ENDDO
	775	C ELSE
	776	C DO ij=1,iip1
	777	C dyqmax(ij)=pente_max*abs(dyqv(ij))
	778	C ENDDO
	779	C ENDIF
	780	C
	781	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	782	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	783	C &THEN
	784	C DO ij=ip1jm+1,ip1jmp1
	785	C dyqmax(ij)=0.
	786	C ENDDO
	787	C ELSE
	788	C DO ij=ip1jm+1,ip1jmp1
	789	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
	790	C ENDDO
	791	C ENDIF
	792	C fin changement 10 07 96
	793	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	794
	795	c calcul des pentes limitees
	796	ijb=ij_begin-iip1
	797	ije=ij_end+iip1
	798	if (pole_nord) ijb=ij_begin+iip1
	799	if (pole_sud) ije=ij_end-iip1
	800
	801	DO ij=ijb,ije
	802	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
	803	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
	804	ELSE
	805	dyq(ij,l)=0.
	806	ENDIF
	807	ENDDO
	808
	809	ENDDO
[764]	810	c$OMP END DO NOWAIT
[630]	811
	812	ijb=ij_begin-iip1
	813	ije=ij_end
	814	if (pole_nord) ijb=ij_begin
	815	if (pole_sud) ije=ij_end-iip1
	816
[764]	817	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	818	DO l=1,llm
	819	DO ij=ijb,ije
	820	IF(masse_adv_v(ij,l).gt.0) THEN
	821	qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
	822	, 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
	823	ELSE
	824	qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
	825	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
	826	ENDIF
	827	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
	828	ENDDO
	829	ENDDO
[764]	830	c$OMP END DO NOWAIT
[630]	831
	832	ijb=ij_begin
	833	ije=ij_end
	834	if (pole_nord) ijb=ij_begin+iip1
	835	if (pole_sud) ije=ij_end-iip1
	836
[764]	837	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	838	DO l=1,llm
	839	DO ij=ijb,ije
	840	newmasse=masse(ij,l)
	841	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
	842
	843	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
	844	& /newmasse
	845	masse(ij,l)=newmasse
	846	ENDDO
	847	c.-. ancienne version
	848	c convpn=SSUM(iim,qbyv(1,l),1)/apoln
	849	c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
	850	if (pole_nord) then
	851	convpn=SSUM(iim,qbyv(1,l),1)
	852	convmpn=ssum(iim,masse_adv_v(1,l),1)
	853	massepn=ssum(iim,masse(1,l),1)
	854	qpn=0.
	855	do ij=1,iim
	856	qpn=qpn+masse(ij,l)*q(ij,l)
	857	enddo
	858	qpn=(qpn+convpn)/(massepn+convmpn)
	859	do ij=1,iip1
	860	q(ij,l)=qpn
	861	enddo
	862	endif
	863
	864	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
	865	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
	866
	867	if (pole_sud) then
	868
	869	convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	870	convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	871	masseps=ssum(iim, masse(ip1jm+1,l),1)
	872	qps=0.
	873	do ij = ip1jm+1,ip1jmp1-1
	874	qps=qps+masse(ij,l)*q(ij,l)
	875	enddo
	876	qps=(qps+convps)/(masseps+convmps)
	877	do ij=ip1jm+1,ip1jmp1
	878	q(ij,l)=qps
	879	enddo
	880	endif
	881	c.-. fin ancienne version
	882
	883	c._. nouvelle version
	884	c convpn=SSUM(iim,qbyv(1,l),1)
	885	c convmpn=ssum(iim,masse_adv_v(1,l),1)
	886	c oldmasse=ssum(iim,masse(1,l),1)
	887	c newmasse=oldmasse+convmpn
	888	c newq=(q(1,l)*oldmasse+convpn)/newmasse
	889	c newmasse=newmasse/apoln
	890	c DO ij = 1,iip1
	891	c q(ij,l)=newq
	892	c masse(ij,l)=newmasse*aire(ij)
	893	c ENDDO
	894	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	895	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	896	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
	897	c newmasse=oldmasse+convmps
	898	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
	899	c newmasse=newmasse/apols
	900	c DO ij = ip1jm+1,ip1jmp1
	901	c q(ij,l)=newq
	902	c masse(ij,l)=newmasse*aire(ij)
	903	c ENDDO
	904	c._. fin nouvelle version
	905	ENDDO
[764]	906	c$OMP END DO NOWAIT
[630]	907
	908	RETURN
	909	END
	910
	911
	912
	913	SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
	914	c
	915	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	916	c
	917	c ********************************************************************
	918	c Shema d'advection " pseudo amont " .
	919	c ********************************************************************
	920	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	921	c dq sont des arguments de sortie pour le s-pg ....
	922	c
	923	c
	924	c --------------------------------------------------------------------
	925	USE Parallel
	926	IMPLICIT NONE
	927	c
	928	#include "dimensions.h"
	929	#include "paramet.h"
	930	#include "logic.h"
	931	#include "comvert.h"
	932	#include "comconst.h"
	933	c
	934	c
	935	c Arguments:
	936	c ----------
	937	REAL masse(ip1jmp1,llm),pente_max
	938	REAL q(ip1jmp1,llm)
	939	REAL w(ip1jmp1,llm+1)
	940	c
	941	c Local
	942	c ---------
	943	c
	944	INTEGER i,ij,l,j,ii
	945	c
[764]	946	REAL,SAVE :: wq(ip1jmp1,llm+1)
	947	REAL newmasse
[630]	948
[764]	949	REAL,SAVE :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm)
	950	REAL dzqmax
[630]	951	REAL sigw
	952
	953	LOGICAL testcpu
	954	SAVE testcpu
[764]	955	c$OMP THREADPRIVATE(testcpu)
[630]	956	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
	957	SAVE temps0,temps1,temps2,temps3,temps4,temps5
[764]	958	c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
	959
[630]	960	REAL SSUM
[764]	961	EXTERNAL SSUM
[630]	962
	963	DATA testcpu/.false./
	964	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
	965	INTEGER ijb,ije,ijb_x,ije_x
	966	c On oriente tout dans le sens de la pression c'est a dire dans le
	967	c sens de W
	968
	969	#ifdef BIDON
	970	IF(testcpu) THEN
	971	temps0=second(0.)
	972	ENDIF
	973	#endif
	974
	975	ijb=ijb_x
	976	ije=ije_x
[764]	977
	978	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	979	DO l=2,llm
	980	DO ij=ijb,ije
	981	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
	982	adzqw(ij,l)=abs(dzqw(ij,l))
	983	ENDDO
	984	ENDDO
[764]	985	c$OMP END DO
[630]	986
[764]	987	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	988	DO l=2,llm-1
	989	DO ij=ijb,ije
	990	#ifdef CRAY
	991	dzq(ij,l)=0.5*
	992	, cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
	993	#else
	994	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
	995	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
	996	ELSE
	997	dzq(ij,l)=0.
	998	ENDIF
	999	#endif
	1000	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
	1001	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
	1002	ENDDO
	1003	ENDDO
[764]	1004	c$OMP END DO NOWAIT
[630]	1005
[764]	1006	c$OMP MASTER
[630]	1007	DO ij=ijb,ije
	1008	dzq(ij,1)=0.
	1009	dzq(ij,llm)=0.
	1010	ENDDO
[764]	1011	c$OMP END MASTER
	1012	c$OMP BARRIER
[630]	1013	#ifdef BIDON
	1014	IF(testcpu) THEN
	1015	temps1=temps1+second(0.)-temps0
	1016	ENDIF
	1017	#endif
	1018	c ---------------------------------------------------------------
	1019	c .... calcul des termes d'advection verticale .......
	1020	c ---------------------------------------------------------------
	1021
	1022	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
	1023
[764]	1024	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1025	DO l = 1,llm-1
	1026	do ij = ijb,ije
	1027	IF(w(ij,l+1).gt.0.) THEN
	1028	sigw=w(ij,l+1)/masse(ij,l+1)
	1029	wq(ij,l+1)=w(ij,l+1)(q(ij,l+1)+0.5(1.-sigw)*dzq(ij,l+1))
	1030	ELSE
	1031	sigw=w(ij,l+1)/masse(ij,l)
	1032	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
	1033	ENDIF
	1034	ENDDO
	1035	ENDDO
[764]	1036	c$OMP END DO NOWAIT
[630]	1037
[764]	1038	c$OMP MASTER
[630]	1039	DO ij=ijb,ije
	1040	wq(ij,llm+1)=0.
	1041	wq(ij,1)=0.
	1042	ENDDO
[764]	1043	c$OMP END MASTER
	1044	c$OMP BARRIER
[630]	1045
[764]	1046	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	1047	DO l=1,llm
	1048	DO ij=ijb,ije
	1049	newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
	1050	q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
	1051	& /newmasse
	1052	masse(ij,l)=newmasse
	1053	ENDDO
	1054	ENDDO
[764]	1055	c$OMP END DO NOWAIT
[630]	1056
	1057
	1058	RETURN
	1059	END
	1060	c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
	1061	c
	1062	c#include "dimensions.h"
	1063	c#include "paramet.h"
	1064
	1065	c CHARACTER() comment
	1066	c real qmin,qmax
	1067	c real zq(ip1jmp1,llm)
	1068
	1069	c INTEGER jadrs(ip1jmp1), jbad, k, i
	1070
	1071
	1072	c DO k = 1, llm
	1073	c jbad = 0
	1074	c DO i = 1, ip1jmp1
	1075	c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
	1076	c jbad = jbad + 1
	1077	c jadrs(jbad) = i
	1078	c ENDIF
	1079	c ENDDO
	1080	c IF (jbad.GT.0) THEN
	1081	c PRINT*, comment
	1082	c DO i = 1, jbad
	1083	cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
	1084	c ENDDO
	1085	c ENDIF
	1086	c ENDDO
	1087
	1088	c return
	1089	c end
	1090
	1091
	1092	subroutine minmaxq_p(zq,qmin,qmax,comment)
	1093
	1094	#include "dimensions.h"
	1095	#include "paramet.h"
	1096
	1097	character*20 comment
	1098	real qmin,qmax
	1099	real zq(ip1jmp1,llm)
	1100	real zzq(iip1,jjp1,llm)
	1101
	1102	integer imin,jmin,lmin,ijlmin
	1103	integer imax,jmax,lmax,ijlmax
	1104
	1105	integer ismin,ismax
	1106
	1107	#ifdef isminismax
	1108	call scopy (ip1jmp1*llm,zq,1,zzq,1)
	1109
	1110	ijlmin=ismin(ijp1llm,zq,1)
	1111	lmin=(ijlmin-1)/ip1jmp1+1
	1112	ijlmin=ijlmin-(lmin-1.)*ip1jmp1
	1113	jmin=(ijlmin-1)/iip1+1
	1114	imin=ijlmin-(jmin-1.)*iip1
	1115	zqmin=zq(ijlmin,lmin)
	1116
	1117	ijlmax=ismax(ijp1llm,zq,1)
	1118	lmax=(ijlmax-1)/ip1jmp1+1
	1119	ijlmax=ijlmax-(lmax-1.)*ip1jmp1
	1120	jmax=(ijlmax-1)/iip1+1
	1121	imax=ijlmax-(jmax-1.)*iip1
	1122	zqmax=zq(ijlmax,lmax)
	1123
	1124	if(zqmin.lt.qmin)
	1125	c s write(*,9999) comment,
	1126	s write(,) comment,
	1127	s imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
	1128	if(zqmax.gt.qmax)
	1129	c s write(*,9999) comment,
	1130	s write(,) comment,
	1131	s imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
	1132	#endif
	1133	return
	1134	9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
	1135	end
	1136
	1137
	1138

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

Download in other formats:

Original Format