Context Navigation

vlsplt_p.F @ 701

Last change on this file since 701 was 630, checked in by Laurent Fairhead, 20 years ago
Import d'une version parallele de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.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
	230	Logical extremum,first,testcpu
	231	SAVE first,testcpu
	232
	233	REAL SSUM
	234	EXTERNAL SSUM
	235	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
	236	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	237
	238	REAL z1,z2,z3
	239
	240	DATA first,testcpu/.true.,.false./
	241	INTEGER ijb,ije,ijb_x,ije_x
	242
	243	IF(first) THEN
	244	temps1=0.
	245	temps2=0.
	246	temps3=0.
	247	temps4=0.
	248	temps5=0.
	249	first=.false.
	250	ENDIF
	251
	252	c calcul de la pente a droite et a gauche de la maille
	253
	254	ijb=ijb_x
	255	ije=ije_x
	256
	257	if (pole_nord.and.ijb==1) ijb=ijb+iip1
	258	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
	259
	260	IF (pente_max.gt.-1.e-5) THEN
	261	c IF (pente_max.gt.10) THEN
	262
	263	c calcul des pentes avec limitation, Van Leer scheme I:
	264	c -----------------------------------------------------
	265
	266	c calcul de la pente aux points u
	267
	268	DO l = 1, llm
	269
	270	DO ij=ijb,ije-1
	271	dxqu(ij)=q(ij+1,l)-q(ij,l)
	272	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	273	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	274	ENDDO
	275	DO ij=ijb+iip1-1,ije,iip1
	276	dxqu(ij)=dxqu(ij-iim)
	277	c sigu(ij)=sigu(ij-iim)
	278	ENDDO
	279
	280	DO ij=ijb,ije
	281	adxqu(ij)=abs(dxqu(ij))
	282	ENDDO
	283
	284	c calcul de la pente maximum dans la maille en valeur absolue
	285
	286	DO ij=ijb+1,ije
	287	dxqmax(ij,l)=pente_max*
	288	, min(adxqu(ij-1),adxqu(ij))
	289	c limitation subtile
	290	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	291
	292
	293	ENDDO
	294
	295	DO ij=ijb+iip1-1,ije,iip1
	296	dxqmax(ij-iim,l)=dxqmax(ij,l)
	297	ENDDO
	298
	299	DO ij=ijb+1,ije
	300	#ifdef CRAY
	301	dxq(ij,l)=
	302	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	303	#else
	304	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	305	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	306	ELSE
	307	c extremum local
	308	dxq(ij,l)=0.
	309	ENDIF
	310	#endif
	311	dxq(ij,l)=0.5*dxq(ij,l)
	312	dxq(ij,l)=
	313	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	314	ENDDO
	315
	316	ENDDO ! l=1,llm
	317	c print*,'Ok calcul des pentes'
	318
	319	ELSE ! (pente_max.lt.-1.e-5)
	320
	321	c Pentes produits:
	322	c ----------------
	323
	324	DO l = 1, llm
	325	DO ij=ijb,ije-1
	326	dxqu(ij)=q(ij+1,l)-q(ij,l)
	327	ENDDO
	328	DO ij=ijb+iip1-1,ije,iip1
	329	dxqu(ij)=dxqu(ij-iim)
	330	ENDDO
	331
	332	DO ij=ijb+1,ije
	333	zz(ij)=dxqu(ij-1)*dxqu(ij)
	334	zz(ij)=zz(ij)+zz(ij)
	335	IF(zz(ij).gt.0) THEN
	336	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	337	ELSE
	338	c extremum local
	339	dxq(ij,l)=0.
	340	ENDIF
	341	ENDDO
	342
	343	ENDDO
	344
	345	ENDIF ! (pente_max.lt.-1.e-5)
	346
	347	c bouclage de la pente en iip1:
	348	c -----------------------------
	349
	350	DO l=1,llm
	351	DO ij=ijb+iip1-1,ije,iip1
	352	dxq(ij-iim,l)=dxq(ij,l)
	353	ENDDO
	354	DO ij=1,ip1jmp1
	355	iadvplus(ij,l)=0
	356	ENDDO
	357
	358	ENDDO
	359
	360	c print*,'Bouclage en iip1'
	361
	362	c calcul des flux a gauche et a droite
	363
	364	#ifdef CRAY
	365
	366	DO l=1,llm
	367	DO ij=ijb,ije-1
	368	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	369	, 1.+u_m(ij,l)/masse(ij+1,l),
	370	, u_m(ij,l))
	371	zdum(ij,l)=0.5*zdum(ij,l)
	372	u_mq(ij,l)=cvmgp(
	373	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	374	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	375	, u_m(ij,l))
	376	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	377	ENDDO
	378	ENDDO
	379	#else
	380	c on cumule le flux correspondant a toutes les mailles dont la masse
	381	c au travers de la paroi pENDant le pas de temps.
	382	c print*,'Cumule ....'
	383
	384	DO l=1,llm
	385	DO ij=ijb,ije-1
	386	c print*,'masse(',ij,')=',masse(ij,l)
	387	IF (u_m(ij,l).gt.0.) THEN
	388	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	389	u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
	390	ELSE
	391	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	392	u_mq(ij,l)=u_m(ij,l)(q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
	393	ENDIF
	394	ENDDO
	395	ENDDO
	396	#endif
	397	c stop
	398
	399	c go to 9999
	400	c detection des points ou on advecte plus que la masse de la
	401	c maille
	402	DO l=1,llm
	403	DO ij=ijb,ije-1
	404	IF(zdum(ij,l).lt.0) THEN
	405	iadvplus(ij,l)=1
	406	u_mq(ij,l)=0.
	407	ENDIF
	408	ENDDO
	409	ENDDO
	410	c print*,'Ok test 1'
	411	DO l=1,llm
	412	DO ij=ijb+iip1-1,ije,iip1
	413	iadvplus(ij,l)=iadvplus(ij-iim,l)
	414	ENDDO
	415	ENDDO
	416	c print*,'Ok test 2'
	417
	418
	419	c traitement special pour le cas ou on advecte en longitude plus que le
	420	c contenu de la maille.
	421	c cette partie est mal vectorisee.
	422
	423	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	424
	425	n0=0
	426	DO l=1,llm
	427	nl(l)=0
	428	DO ij=ijb,ije
	429	nl(l)=nl(l)+iadvplus(ij,l)
	430	ENDDO
	431	n0=n0+nl(l)
	432	ENDDO
	433
	434	cym IF(n0.gt.1) THEN
	435	IF(n0.gt.0) THEN
	436
	437	c PRINT*,'Nombre de points pour lesquels on advect plus que le'
	438	c & ,'contenu de la maille : ',n0
	439
	440	DO l=1,llm
	441	IF(nl(l).gt.0) THEN
	442	iju=0
	443	c indicage des mailles concernees par le traitement special
	444	DO ij=ijb,ije
	445	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	446	iju=iju+1
	447	indu(iju)=ij
	448	ENDIF
	449	ENDDO
	450	niju=iju
	451	c PRINT*,'niju,nl',niju,nl(l)
	452
	453	c traitement des mailles
	454	DO iju=1,niju
	455	ij=indu(iju)
	456	j=(ij-1)/iip1+1
	457	zu_m=u_m(ij,l)
	458	u_mq(ij,l)=0.
	459	IF(zu_m.gt.0.) THEN
	460	ijq=ij
	461	i=ijq-(j-1)*iip1
	462	c accumulation pour les mailles completements advectees
	463	do while(zu_m.gt.masse(ijq,l))
	464	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	465	zu_m=zu_m-masse(ijq,l)
	466	i=mod(i-2+iim,iim)+1
	467	ijq=(j-1)*iip1+i
	468	ENDDO
	469	c ajout de la maille non completement advectee
	470	u_mq(ij,l)=u_mq(ij,l)+zu_m*
	471	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
	472	ELSE
	473	ijq=ij+1
	474	i=ijq-(j-1)*iip1
	475	c accumulation pour les mailles completements advectees
	476	do while(-zu_m.gt.masse(ijq,l))
	477	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	478	zu_m=zu_m+masse(ijq,l)
	479	i=mod(i,iim)+1
	480	ijq=(j-1)*iip1+i
	481	ENDDO
	482	c ajout de la maille non completement advectee
	483	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
	484	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
	485	ENDIF
	486	ENDDO
	487	ENDIF
	488	ENDDO
	489	ENDIF ! n0.gt.0
	490	9999 continue
	491
	492
	493	c bouclage en latitude
	494	c print*,'Avant bouclage en latitude'
	495	DO l=1,llm
	496	DO ij=ijb+iip1-1,ije,iip1
	497	u_mq(ij,l)=u_mq(ij-iim,l)
	498	ENDDO
	499	ENDDO
	500
	501
	502	c calcul des tENDances
	503
	504	DO l=1,llm
	505	DO ij=ijb+1,ije
	506	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	507	q(ij,l)=(q(ij,l)*masse(ij,l)+
	508	& u_mq(ij-1,l)-u_mq(ij,l))
	509	& /new_m
	510	masse(ij,l)=new_m
	511	ENDDO
	512	c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	513	DO ij=ijb+iip1-1,ije,iip1
	514	q(ij-iim,l)=q(ij,l)
	515	masse(ij-iim,l)=masse(ij,l)
	516	ENDDO
	517	ENDDO
	518	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	519	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	520
	521
	522	RETURN
	523	END
	524
	525
	526	SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
	527	c
	528	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	529	c
	530	c ********************************************************************
	531	c Shema d'advection " pseudo amont " .
	532	c ********************************************************************
	533	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	534	c dq sont des arguments de sortie pour le s-pg ....
	535	c
	536	c
	537	c --------------------------------------------------------------------
	538	USE parallel
	539	IMPLICIT NONE
	540	c
	541	#include "dimensions.h"
	542	#include "paramet.h"
	543	#include "logic.h"
	544	#include "comvert.h"
	545	#include "comconst.h"
	546	#include "comgeom.h"
	547	c
	548	c
	549	c Arguments:
	550	c ----------
	551	REAL masse(ip1jmp1,llm),pente_max
	552	REAL masse_adv_v( ip1jm,llm)
	553	REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
	554	c
	555	c Local
	556	c ---------
	557	c
	558	INTEGER i,ij,l
	559	c
	560	REAL airej2,airejjm,airescb(iim),airesch(iim)
	561	REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
	562	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
	563	REAL qbyv(ip1jm,llm)
	564
	565	REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	566	c REAL newq,oldmasse
	567	Logical extremum,first,testcpu
	568	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
	569	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	570	SAVE first,testcpu
	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
	577	SAVE airej2,airejjm
	578	c
	579	c
	580	REAL SSUM
	581	EXTERNAL SSUM
	582
	583	DATA first,testcpu/.true.,.false./
	584	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
	585	INTEGER ijb,ije
	586
	587	IF(first) THEN
	588	c PRINT*,'Shema Amont nouveau appele dans Vanleer '
	589	first=.false.
	590	do i=2,iip1
	591	coslon(i)=cos(rlonv(i))
	592	sinlon(i)=sin(rlonv(i))
	593	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	594	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	595	ENDDO
	596	coslon(1)=coslon(iip1)
	597	coslondlon(1)=coslondlon(iip1)
	598	sinlon(1)=sinlon(iip1)
	599	sinlondlon(1)=sinlondlon(iip1)
	600	airej2 = SSUM( iim, aire(iip2), 1 )
	601	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
	602	ENDIF
	603
	604	c
	605	c PRINT*,'CALCUL EN LATITUDE'
	606
	607
	608	DO l = 1, llm
	609	c
	610	c --------------------------------
	611	c CALCUL EN LATITUDE
	612	c --------------------------------
	613
	614	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
	615	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
	616	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
	617
	618	if (pole_nord) then
	619	DO i = 1, iim
	620	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
	621	ENDDO
	622	qpns = SSUM( iim, airescb ,1 ) / airej2
	623	endif
	624
	625	if (pole_sud) then
	626	DO i = 1, iim
	627	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
	628	ENDDO
	629	qpsn = SSUM( iim, airesch ,1 ) / airejjm
	630	endif
	631
	632
	633
	634	c calcul des pentes aux points v
	635
	636	ijb=ij_begin-2*iip1
	637	ije=ij_end+iip1
	638	if (pole_nord) ijb=ij_begin
	639	if (pole_sud) ije=ij_end-iip1
	640
	641	DO ij=ijb,ije
	642	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
	643	adyqv(ij)=abs(dyqv(ij))
	644	ENDDO
	645
	646	c calcul des pentes aux points scalaires
	647	ijb=ij_begin-iip1
	648	ije=ij_end+iip1
	649	if (pole_nord) ijb=ij_begin+iip1
	650	if (pole_sud) ije=ij_end-iip1
	651
	652	DO ij=ijb,ije
	653	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
	654	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
	655	dyqmax(ij)=pente_max*dyqmax(ij)
	656	ENDDO
	657
	658	c calcul des pentes aux poles
	659	IF (pole_nord) THEN
	660	DO ij=1,iip1
	661	dyq(ij,l)=qpns-q(ij+iip1,l)
	662	ENDDO
	663
	664	dyn1=0.
	665	dyn2=0.
	666	DO ij=1,iim
	667	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
	668	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
	669	ENDDO
	670	DO ij=1,iip1
	671	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
	672	ENDDO
	673
	674	DO ij=1,iip1
	675	dyq(ij,l)=0.
	676	ENDDO
	677	c ym tout cela ne sert pas a grand chose
	678	ENDIF
	679
	680	IF (pole_sud) THEN
	681
	682	DO ij=1,iip1
	683	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
	684	ENDDO
	685
	686	dys1=0.
	687	dys2=0.
	688
	689	DO ij=1,iim
	690	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
	691	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
	692	ENDDO
	693
	694	DO ij=1,iip1
	695	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
	696	ENDDO
	697
	698	DO ij=1,iip1
	699	dyq(ip1jm+ij,l)=0.
	700	ENDDO
	701	c ym tout cela ne sert pas a grand chose
	702	ENDIF
	703
	704	c filtrage de la derivee
	705
	706	c calcul des pentes limites aux poles
	707	c ym partie inutile
	708	c goto 8888
	709	c fn=1.
	710	c fs=1.
	711	c DO ij=1,iim
	712	c IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
	713	c fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
	714	c ENDIF
	715	c IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
	716	c fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
	717	c ENDIF
	718	c ENDDO
	719	c DO ij=1,iip1
	720	c dyq(ij,l)=fn*dyq(ij,l)
	721	c dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
	722	c ENDDO
	723	c 8888 continue
	724
	725
	726	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	727	C En memoire de dIFferents tests sur la
	728	C limitation des pentes aux poles.
	729	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	730	C PRINT*,dyq(1)
	731	C PRINT*,dyqv(iip1+1)
	732	C apn=abs(dyq(1)/dyqv(iip1+1))
	733	C PRINT*,dyq(ip1jm+1)
	734	C PRINT*,dyqv(ip1jm-iip1+1)
	735	C aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
	736	C DO ij=2,iim
	737	C apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
	738	C aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
	739	C ENDDO
	740	C apn=min(pente_max/apn,1.)
	741	C aps=min(pente_max/aps,1.)
	742	C
	743	C
	744	C cas ou on a un extremum au pole
	745	C
	746	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	747	C & apn=0.
	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 & aps=0.
	751	C
	752	C limitation des pentes aux poles
	753	C DO ij=1,iip1
	754	C dyq(ij)=apn*dyq(ij)
	755	C dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
	756	C ENDDO
	757	C
	758	C test
	759	C DO ij=1,iip1
	760	C dyq(iip1+ij)=0.
	761	C dyq(ip1jm+ij-iip1)=0.
	762	C ENDDO
	763	C DO ij=1,ip1jmp1
	764	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
	765	C ENDDO
	766	C
	767	C changement 10 07 96
	768	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	769	C & THEN
	770	C DO ij=1,iip1
	771	C dyqmax(ij)=0.
	772	C ENDDO
	773	C ELSE
	774	C DO ij=1,iip1
	775	C dyqmax(ij)=pente_max*abs(dyqv(ij))
	776	C ENDDO
	777	C ENDIF
	778	C
	779	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	780	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	781	C &THEN
	782	C DO ij=ip1jm+1,ip1jmp1
	783	C dyqmax(ij)=0.
	784	C ENDDO
	785	C ELSE
	786	C DO ij=ip1jm+1,ip1jmp1
	787	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
	788	C ENDDO
	789	C ENDIF
	790	C fin changement 10 07 96
	791	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	792
	793	c calcul des pentes limitees
	794	ijb=ij_begin-iip1
	795	ije=ij_end+iip1
	796	if (pole_nord) ijb=ij_begin+iip1
	797	if (pole_sud) ije=ij_end-iip1
	798
	799	DO ij=ijb,ije
	800	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
	801	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
	802	ELSE
	803	dyq(ij,l)=0.
	804	ENDIF
	805	ENDDO
	806
	807	ENDDO
	808
	809	ijb=ij_begin-iip1
	810	ije=ij_end
	811	if (pole_nord) ijb=ij_begin
	812	if (pole_sud) ije=ij_end-iip1
	813
	814
	815	DO l=1,llm
	816	DO ij=ijb,ije
	817	IF(masse_adv_v(ij,l).gt.0) THEN
	818	qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
	819	, 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
	820	ELSE
	821	qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
	822	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
	823	ENDIF
	824	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
	825	ENDDO
	826	ENDDO
	827
	828
	829	ijb=ij_begin
	830	ije=ij_end
	831	if (pole_nord) ijb=ij_begin+iip1
	832	if (pole_sud) ije=ij_end-iip1
	833
	834	DO l=1,llm
	835	DO ij=ijb,ije
	836	newmasse=masse(ij,l)
	837	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
	838
	839	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
	840	& /newmasse
	841	masse(ij,l)=newmasse
	842	ENDDO
	843	c.-. ancienne version
	844	c convpn=SSUM(iim,qbyv(1,l),1)/apoln
	845	c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
	846	if (pole_nord) then
	847	convpn=SSUM(iim,qbyv(1,l),1)
	848	convmpn=ssum(iim,masse_adv_v(1,l),1)
	849	massepn=ssum(iim,masse(1,l),1)
	850	qpn=0.
	851	do ij=1,iim
	852	qpn=qpn+masse(ij,l)*q(ij,l)
	853	enddo
	854	qpn=(qpn+convpn)/(massepn+convmpn)
	855	do ij=1,iip1
	856	q(ij,l)=qpn
	857	enddo
	858	endif
	859
	860	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
	861	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
	862
	863	if (pole_sud) then
	864
	865	convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	866	convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	867	masseps=ssum(iim, masse(ip1jm+1,l),1)
	868	qps=0.
	869	do ij = ip1jm+1,ip1jmp1-1
	870	qps=qps+masse(ij,l)*q(ij,l)
	871	enddo
	872	qps=(qps+convps)/(masseps+convmps)
	873	do ij=ip1jm+1,ip1jmp1
	874	q(ij,l)=qps
	875	enddo
	876	endif
	877	c.-. fin ancienne version
	878
	879	c._. nouvelle version
	880	c convpn=SSUM(iim,qbyv(1,l),1)
	881	c convmpn=ssum(iim,masse_adv_v(1,l),1)
	882	c oldmasse=ssum(iim,masse(1,l),1)
	883	c newmasse=oldmasse+convmpn
	884	c newq=(q(1,l)*oldmasse+convpn)/newmasse
	885	c newmasse=newmasse/apoln
	886	c DO ij = 1,iip1
	887	c q(ij,l)=newq
	888	c masse(ij,l)=newmasse*aire(ij)
	889	c ENDDO
	890	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	891	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	892	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
	893	c newmasse=oldmasse+convmps
	894	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
	895	c newmasse=newmasse/apols
	896	c DO ij = ip1jm+1,ip1jmp1
	897	c q(ij,l)=newq
	898	c masse(ij,l)=newmasse*aire(ij)
	899	c ENDDO
	900	c._. fin nouvelle version
	901	ENDDO
	902
	903	RETURN
	904	END
	905
	906
	907
	908	SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
	909	c
	910	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	911	c
	912	c ********************************************************************
	913	c Shema d'advection " pseudo amont " .
	914	c ********************************************************************
	915	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	916	c dq sont des arguments de sortie pour le s-pg ....
	917	c
	918	c
	919	c --------------------------------------------------------------------
	920	USE Parallel
	921	IMPLICIT NONE
	922	c
	923	#include "dimensions.h"
	924	#include "paramet.h"
	925	#include "logic.h"
	926	#include "comvert.h"
	927	#include "comconst.h"
	928	c
	929	c
	930	c Arguments:
	931	c ----------
	932	REAL masse(ip1jmp1,llm),pente_max
	933	REAL q(ip1jmp1,llm)
	934	REAL w(ip1jmp1,llm+1)
	935	c
	936	c Local
	937	c ---------
	938	c
	939	INTEGER i,ij,l,j,ii
	940	c
	941	REAL wq(ip1jmp1,llm+1),newmasse
	942
	943	REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
	944	REAL sigw
	945
	946	LOGICAL testcpu
	947	SAVE testcpu
	948
	949	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
	950	SAVE temps0,temps1,temps2,temps3,temps4,temps5
	951	REAL SSUM
	952	EXTERNAL SSUM, convflu
	953	EXTERNAL filtreg
	954
	955	DATA testcpu/.false./
	956	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
	957	INTEGER ijb,ije,ijb_x,ije_x
	958	c On oriente tout dans le sens de la pression c'est a dire dans le
	959	c sens de W
	960
	961	#ifdef BIDON
	962	IF(testcpu) THEN
	963	temps0=second(0.)
	964	ENDIF
	965	#endif
	966
	967	ijb=ijb_x
	968	ije=ije_x
	969
	970	DO l=2,llm
	971	DO ij=ijb,ije
	972	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
	973	adzqw(ij,l)=abs(dzqw(ij,l))
	974	ENDDO
	975	ENDDO
	976
	977	DO l=2,llm-1
	978	DO ij=ijb,ije
	979	#ifdef CRAY
	980	dzq(ij,l)=0.5*
	981	, cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
	982	#else
	983	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
	984	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
	985	ELSE
	986	dzq(ij,l)=0.
	987	ENDIF
	988	#endif
	989	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
	990	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
	991	ENDDO
	992	ENDDO
	993
	994	DO ij=ijb,ije
	995	dzq(ij,1)=0.
	996	dzq(ij,llm)=0.
	997	ENDDO
	998
	999	#ifdef BIDON
	1000	IF(testcpu) THEN
	1001	temps1=temps1+second(0.)-temps0
	1002	ENDIF
	1003	#endif
	1004	c ---------------------------------------------------------------
	1005	c .... calcul des termes d'advection verticale .......
	1006	c ---------------------------------------------------------------
	1007
	1008	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
	1009
	1010	DO l = 1,llm-1
	1011	do ij = ijb,ije
	1012	IF(w(ij,l+1).gt.0.) THEN
	1013	sigw=w(ij,l+1)/masse(ij,l+1)
	1014	wq(ij,l+1)=w(ij,l+1)(q(ij,l+1)+0.5(1.-sigw)*dzq(ij,l+1))
	1015	ELSE
	1016	sigw=w(ij,l+1)/masse(ij,l)
	1017	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
	1018	ENDIF
	1019	ENDDO
	1020	ENDDO
	1021
	1022	DO ij=ijb,ije
	1023	wq(ij,llm+1)=0.
	1024	wq(ij,1)=0.
	1025	ENDDO
	1026
	1027	DO l=1,llm
	1028	DO ij=ijb,ije
	1029	newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
	1030	q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
	1031	& /newmasse
	1032	masse(ij,l)=newmasse
	1033	ENDDO
	1034	ENDDO
	1035
	1036
	1037	RETURN
	1038	END
	1039	c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
	1040	c
	1041	c#include "dimensions.h"
	1042	c#include "paramet.h"
	1043
	1044	c CHARACTER() comment
	1045	c real qmin,qmax
	1046	c real zq(ip1jmp1,llm)
	1047
	1048	c INTEGER jadrs(ip1jmp1), jbad, k, i
	1049
	1050
	1051	c DO k = 1, llm
	1052	c jbad = 0
	1053	c DO i = 1, ip1jmp1
	1054	c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
	1055	c jbad = jbad + 1
	1056	c jadrs(jbad) = i
	1057	c ENDIF
	1058	c ENDDO
	1059	c IF (jbad.GT.0) THEN
	1060	c PRINT*, comment
	1061	c DO i = 1, jbad
	1062	cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
	1063	c ENDDO
	1064	c ENDIF
	1065	c ENDDO
	1066
	1067	c return
	1068	c end
	1069
	1070
	1071	subroutine minmaxq_p(zq,qmin,qmax,comment)
	1072
	1073	#include "dimensions.h"
	1074	#include "paramet.h"
	1075
	1076	character*20 comment
	1077	real qmin,qmax
	1078	real zq(ip1jmp1,llm)
	1079	real zzq(iip1,jjp1,llm)
	1080
	1081	integer imin,jmin,lmin,ijlmin
	1082	integer imax,jmax,lmax,ijlmax
	1083
	1084	integer ismin,ismax
	1085
	1086	#ifdef isminismax
	1087	call scopy (ip1jmp1*llm,zq,1,zzq,1)
	1088
	1089	ijlmin=ismin(ijp1llm,zq,1)
	1090	lmin=(ijlmin-1)/ip1jmp1+1
	1091	ijlmin=ijlmin-(lmin-1.)*ip1jmp1
	1092	jmin=(ijlmin-1)/iip1+1
	1093	imin=ijlmin-(jmin-1.)*iip1
	1094	zqmin=zq(ijlmin,lmin)
	1095
	1096	ijlmax=ismax(ijp1llm,zq,1)
	1097	lmax=(ijlmax-1)/ip1jmp1+1
	1098	ijlmax=ijlmax-(lmax-1.)*ip1jmp1
	1099	jmax=(ijlmax-1)/iip1+1
	1100	imax=ijlmax-(jmax-1.)*iip1
	1101	zqmax=zq(ijlmax,lmax)
	1102
	1103	if(zqmin.lt.qmin)
	1104	c s write(*,9999) comment,
	1105	s write(,) comment,
	1106	s imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
	1107	if(zqmax.gt.qmax)
	1108	c s write(*,9999) comment,
	1109	s write(,) comment,
	1110	s imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
	1111	#endif
	1112	return
	1113	9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
	1114	end
	1115
	1116
	1117

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

Download in other formats:

Original Format