Context Navigation

advn.F @ 5353

Last change on this file since 5353 was 774, checked in by Laurent Fairhead, 18 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.9 KB

Rev	Line
[630]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
	5	c
	6	c Auteur : F. Hourdin
	7	c
	8	c ********************************************************************
	9	c Shema d'advection " pseudo amont " .
	10	c ********************************************************************
	11	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	12	c
	13	c pbaru,pbarv,w flux de masse en u ,v ,w
	14	c pdt pas de temps
	15	c
	16	c --------------------------------------------------------------------
	17	IMPLICIT NONE
	18	c
	19	#include "dimensions.h"
	20	#include "paramet.h"
	21	#include "logic.h"
	22	#include "comvert.h"
	23	#include "comconst.h"
	24	#include "comgeom.h"
	25	#include "iniprint.h"
	26
	27	c
	28	c Arguments:
	29	c ----------
	30	integer mode
	31	real masse(ip1jmp1,llm)
	32	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
	33	REAL q(ip1jmp1,llm)
	34	REAL w(ip1jmp1,llm),pdt
	35	c
	36	c Local
	37	c ---------
	38	c
	39	INTEGER i,ij,l,j,ii
	40	integer ijlqmin,iqmin,jqmin,lqmin
	41	integer ismin
	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,qpn,qps
	48	real zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)
	49	real zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)
	50	real zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)
	51	real temps0,temps1,temps2,temps3
	52	real ztemps1,ztemps2,ztemps3,ssum
	53	logical testcpu
	54	save testcpu
	55	save temps1,temps2,temps3
	56	real zzpbar,zzw
	57
	58	#ifdef CRAY
	59	real second
	60	#endif
	61
	62	real qmin,qmax
	63	data qmin,qmax/0.,1./
	64	data testcpu/.false./
	65	data temps1,temps2,temps3/0.,0.,0./
	66
	67	zzpbar = 0.5 * pdt
	68	zzw = pdt
	69
	70	DO l=1,llm
	71	DO ij = iip2,ip1jm
	72	mu(ij,l)=pbaru(ij,l) * zzpbar
	73	ENDDO
	74	DO ij=1,ip1jm
	75	mv(ij,l)=pbarv(ij,l) * zzpbar
	76	ENDDO
	77	DO ij=1,ip1jmp1
	78	mw(ij,l)=w(ij,l) * zzw
	79	ENDDO
	80	ENDDO
	81
	82	DO ij=1,ip1jmp1
	83	mw(ij,llm+1)=0.
	84	ENDDO
	85
	86	do l=1,llm
	87	qpn=0.
	88	qps=0.
	89	do ij=1,iim
	90	qpn=qpn+q(ij,l)*masse(ij,l)
	91	qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)
	92	enddo
	93	qpn=qpn/ssum(iim,masse(1,l),1)
	94	qps=qps/ssum(iim,masse(ip1jm+1,l),1)
	95	do ij=1,iip1
	96	q(ij,l)=qpn
	97	q(ip1jm+ij,l)=qps
	98	enddo
	99	enddo
	100
	101	do ij=1,ip1jmp1
	102	mw(ij,llm+1)=0.
	103	enddo
	104	do l=1,llm
	105	do ij=1,ip1jmp1
	106	zq(ij,l)=q(ij,l)
	107	zm(ij,l)=masse(ij,l)
	108	enddo
	109	enddo
	110
	111	c call minmaxq(zq,qmin,qmax,'avant vlx ')
	112	call advnqx(zq,zqg,zqd)
	113	call advnx(zq,zqg,zqd,zm,mu,mode)
	114	call advnqy(zq,zqs,zqn)
	115	call advny(zq,zqs,zqn,zm,mv)
	116	call advnqz(zq,zqh,zqb)
	117	call advnz(zq,zqh,zqb,zm,mw)
	118	c call vlz(zq,0.,zm,mw)
	119	call advnqy(zq,zqs,zqn)
	120	call advny(zq,zqs,zqn,zm,mv)
	121	call advnqx(zq,zqg,zqd)
	122	call advnx(zq,zqg,zqd,zm,mu,mode)
	123	c call minmaxq(zq,qmin,qmax,'apres vlx ')
	124
	125	#ifdef CRAY
	126	if(testcpu) then
	127	ztemps1=second(0.)
	128	temps1=temps1+ztemps1-ztemps2
	129	print*,'VLSPLT X:',temps1,' Y:',temps2,' Z:',temps3
	130	endif
	131	#endif
	132	do l=1,llm
	133	do ij=1,ip1jmp1
	134	q(ij,l)=zq(ij,l)
	135	enddo
	136	do ij=1,ip1jm+1,iip1
	137	q(ij+iim,l)=q(ij,l)
	138	enddo
	139	enddo
	140
	141	RETURN
	142	END
	143
	144	SUBROUTINE advnqx(q,qg,qd)
	145	c
	146	c Auteurs: Calcul des valeurs de q aux point u.
	147	c
	148	c --------------------------------------------------------------------
	149	IMPLICIT NONE
	150	c
	151	#include "dimensions.h"
	152	#include "paramet.h"
	153	#include "iniprint.h"
	154	c
	155	c
	156	c Arguments:
	157	c ----------
	158	real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
	159	c
	160	c Local
	161	c ---------
	162	c
	163	INTEGER ij,l
	164	c
	165	real dxqu(ip1jmp1),zqu(ip1jmp1)
	166	real zqmax(ip1jmp1),zqmin(ip1jmp1)
	167	logical extremum(ip1jmp1)
	168
	169	integer mode
	170	save mode
	171	data mode/1/
	172
	173	c calcul des pentes en u:
	174	c -----------------------
	175	if (mode.eq.0) then
	176	do l=1,llm
	177	do ij=1,ip1jm
	178	qd(ij,l)=q(ij,l)
	179	qg(ij,l)=q(ij,l)
	180	enddo
	181	enddo
	182	else
	183	do l = 1, llm
	184	do ij=iip2,ip1jm-1
	185	dxqu(ij)=q(ij+1,l)-q(ij,l)
	186	zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))
	187	enddo
	188	do ij=iip1+iip1,ip1jm,iip1
	189	dxqu(ij)=dxqu(ij-iim)
	190	zqu(ij)=zqu(ij-iim)
	191	enddo
	192	do ij=iip2,ip1jm-1
	193	zqu(ij)=zqu(ij)-dxqu(ij+1)/12.
	194	enddo
	195	do ij=iip1+iip1,ip1jm,iip1
	196	zqu(ij)=zqu(ij-iim)
	197	enddo
	198	do ij=iip2+1,ip1jm
	199	zqu(ij)=zqu(ij)+dxqu(ij-1)/12.
	200	enddo
	201	do ij=iip1+iip1,ip1jm,iip1
	202	zqu(ij-iim)=zqu(ij)
	203	enddo
	204
	205	c calcul des valeurs max et min acceptees aux interfaces
	206
	207	do ij=iip2,ip1jm-1
	208	zqmax(ij)=max(q(ij+1,l),q(ij,l))
	209	zqmin(ij)=min(q(ij+1,l),q(ij,l))
	210	enddo
	211	do ij=iip1+iip1,ip1jm,iip1
	212	zqmax(ij)=zqmax(ij-iim)
	213	zqmin(ij)=zqmin(ij-iim)
	214	enddo
	215	do ij=iip2+1,ip1jm
	216	extremum(ij)=dxqu(ij)*dxqu(ij-1).le.0.
	217	enddo
	218	do ij=iip1+iip1,ip1jm,iip1
	219	extremum(ij-iim)=extremum(ij)
	220	enddo
	221	do ij=iip2,ip1jm
	222	zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))
	223	enddo
	224	do ij=iip2+1,ip1jm
	225	if(extremum(ij)) then
	226	qg(ij,l)=q(ij,l)
	227	qd(ij,l)=q(ij,l)
	228	else
	229	qd(ij,l)=zqu(ij)
	230	qg(ij,l)=zqu(ij-1)
	231	endif
	232	enddo
	233	do ij=iip1+iip1,ip1jm,iip1
	234	qd(ij-iim,l)=qd(ij,l)
	235	qg(ij-iim,l)=qg(ij,l)
	236	enddo
	237
	238	goto 8888
	239
	240	do ij=iip2+1,ip1jm
	241	if(extremum(ij).and..not.extremum(ij-1))
	242	s qd(ij-1,l)=q(ij,l)
	243	enddo
	244
	245	do ij=iip1+iip1,ip1jm,iip1
	246	qd(ij-iim,l)=qd(ij,l)
	247	enddo
	248	do ij=iip2,ip1jm-1
	249	if (extremum(ij).and..not.extremum(ij+1))
	250	s qg(ij+1,l)=q(ij,l)
	251	enddo
	252
	253	do ij=iip1+iip1,ip1jm,iip1
	254	qg(ij,l)=qg(ij-iim,l)
	255	enddo
	256	8888 continue
	257	enddo
	258	endif
	259	RETURN
	260	END
	261	SUBROUTINE advnqy(q,qs,qn)
	262	c
	263	c Auteurs: Calcul des valeurs de q aux point v.
	264	c
	265	c --------------------------------------------------------------------
	266	IMPLICIT NONE
	267	c
	268	#include "dimensions.h"
	269	#include "paramet.h"
	270	#include "iniprint.h"
	271	c
	272	c
	273	c Arguments:
	274	c ----------
	275	real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
	276	c
	277	c Local
	278	c ---------
	279	c
	280	INTEGER ij,l
	281	c
	282	real dyqv(ip1jm),zqv(ip1jm,llm)
	283	real zqmax(ip1jm),zqmin(ip1jm)
	284	logical extremum(ip1jmp1)
	285
	286	integer mode
	287	save mode
	288	data mode/1/
	289
	290	if (mode.eq.0) then
	291	do l=1,llm
	292	do ij=1,ip1jmp1
	293	qn(ij,l)=q(ij,l)
	294	qs(ij,l)=q(ij,l)
	295	enddo
	296	enddo
	297	else
	298
	299	c calcul des pentes en u:
	300	c -----------------------
	301	do l = 1, llm
	302	do ij=1,ip1jm
	303	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
	304	enddo
	305
	306	do ij=iip2,ip1jm-iip1
	307	zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))
	308	zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.
	309	enddo
	310
	311	do ij=iip2,ip1jm
	312	extremum(ij)=dyqv(ij)*dyqv(ij-iip1).le.0.
	313	enddo
	314
	315	c Pas de pentes aux poles
	316	do ij=1,iip1
	317	zqv(ij,l)=q(ij,l)
	318	zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)
	319	extremum(ij)=.true.
	320	extremum(ip1jmp1-iip1+ij)=.true.
	321	enddo
	322
	323	c calcul des valeurs max et min acceptees aux interfaces
	324	do ij=1,ip1jm
	325	zqmax(ij)=max(q(ij+iip1,l),q(ij,l))
	326	zqmin(ij)=min(q(ij+iip1,l),q(ij,l))
	327	enddo
	328
	329	do ij=1,ip1jm
	330	zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))
	331	enddo
	332
	333	do ij=iip2,ip1jm
	334	if(extremum(ij)) then
	335	qs(ij,l)=q(ij,l)
	336	qn(ij,l)=q(ij,l)
	337	c if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
	338	c if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
	339	else
	340	qs(ij,l)=zqv(ij,l)
	341	qn(ij,l)=zqv(ij-iip1,l)
	342	endif
	343	enddo
	344
	345	do ij=1,iip1
	346	qs(ij,l)=q(ij,l)
	347	qn(ij,l)=q(ij,l)
	348	qs(ip1jm+ij,l)=q(ip1jm+ij,l)
	349	qn(ip1jm+ij,l)=q(ip1jm+ij,l)
	350	enddo
	351
	352	enddo
	353	endif
	354	RETURN
	355	END
	356
	357	SUBROUTINE advnqz(q,qh,qb)
	358	c
	359	c Auteurs: Calcul des valeurs de q aux point v.
	360	c
	361	c --------------------------------------------------------------------
	362	IMPLICIT NONE
	363	c
	364	#include "dimensions.h"
	365	#include "paramet.h"
	366	#include "iniprint.h"
	367	c
	368	c
	369	c Arguments:
	370	c ----------
	371	real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
	372	c
	373	c Local
	374	c ---------
	375	c
	376	INTEGER ij,l
	377	c
	378	real dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
	379	real zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
	380	logical extremum(ip1jmp1,llm)
	381
	382	integer mode
	383	save mode
	384
	385	data mode/1/
	386
	387	c calcul des pentes en u:
	388	c -----------------------
	389
	390	if (mode.eq.0) then
	391	do l=1,llm
	392	do ij=1,ip1jmp1
	393	qb(ij,l)=q(ij,l)
	394	qh(ij,l)=q(ij,l)
	395	enddo
	396	enddo
	397	else
	398	do l = 2, llm
	399	do ij=1,ip1jmp1
	400	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
	401	zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))
	402	enddo
	403	enddo
	404	do ij=1,ip1jmp1
	405	dzqw(ij,1)=0.
	406	dzqw(ij,llm+1)=0.
	407	enddo
	408	do l=2,llm
	409	do ij=1,ip1jmp1
	410	zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.
	411	enddo
	412	enddo
	413	do l=2,llm-1
	414	do ij=1,ip1jmp1
	415	extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1).le.0.
	416	enddo
	417	enddo
	418
	419	c Pas de pentes en bas et en haut
	420	do ij=1,ip1jmp1
	421	zqw(ij,2)=q(ij,1)
	422	zqw(ij,llm)=q(ij,llm)
	423	extremum(ij,1)=.true.
	424	extremum(ij,llm)=.true.
	425	enddo
	426
	427	c calcul des valeurs max et min acceptees aux interfaces
	428	do l=2,llm
	429	do ij=1,ip1jmp1
	430	zqmax(ij,l)=max(q(ij,l-1),q(ij,l))
	431	zqmin(ij,l)=min(q(ij,l-1),q(ij,l))
	432	enddo
	433	enddo
	434
	435	do l=2,llm
	436	do ij=1,ip1jmp1
	437	zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))
	438	enddo
	439	enddo
	440
	441	do l=2,llm-1
	442	do ij=1,ip1jmp1
	443	if(extremum(ij,l)) then
	444	qh(ij,l)=q(ij,l)
	445	qb(ij,l)=q(ij,l)
	446	else
	447	qh(ij,l)=zqw(ij,l+1)
	448	qb(ij,l)=zqw(ij,l)
	449	endif
	450	enddo
	451	enddo
	452	c do l=2,llm-1
	453	c do ij=1,ip1jmp1
	454	c if(extremum(ij,l)) then
	455	c if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
	456	c if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
	457	c endif
	458	c enddo
	459	c enddo
	460
	461	do ij=1,ip1jmp1
	462	qb(ij,1)=q(ij,1)
	463	qh(ij,1)=q(ij,1)
	464	qb(ij,llm)=q(ij,llm)
	465	qh(ij,llm)=q(ij,llm)
	466	enddo
	467
	468	endif
	469
	470	RETURN
	471	END
	472
	473	SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
	474	c
	475	c Auteur : F. Hourdin
	476	c
	477	c ********************************************************************
	478	c Shema d'advection " pseudo amont " .
	479	c ********************************************************************
	480	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	481	c
	482	c
	483	c --------------------------------------------------------------------
	484	IMPLICIT NONE
	485	c
	486	#include "dimensions.h"
	487	#include "paramet.h"
	488	#include "logic.h"
	489	#include "comvert.h"
	490	#include "comconst.h"
	491	#include "iniprint.h"
	492	c
	493	c
	494	c Arguments:
	495	c ----------
	496	integer mode
	497	real masse(ip1jmp1,llm)
	498	real u_m( ip1jmp1,llm )
	499	real q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
	500	c
	501	c Local
	502	c ---------
	503	c
	504	INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
	505	integer n0,nl(llm)
	506	c
	507	real new_m,zu_m,zdq,zz
	508	real zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig
	509	real u_mq(ip1jmp1,llm)
	510
	511	real zm,zq,zsigm,zsigp,zqm,zqp,zu
	512
	513	logical ladvplus(ip1jmp1,llm)
	514
	515	real prec
	516	save prec
	517
	518	#ifdef CRAY
	519	data prec/1.e-24/
	520	#else
	521	data prec/1.e-15/
	522	#endif
	523
	524	do l=1,llm
	525	do ij=iip2,ip1jm
	526	zdq=qd(ij,l)-qg(ij,l)
	527	c if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
	528	c print*,'probleme au point ij=',ij,' l=',l
	529	c print*,qd(ij,l),q(ij,l),qg(ij,l)
	530	c qd(ij,l)=q(ij,l)
	531	c qg(ij,l)=q(ij,l)
	532	c endif
	533	if(abs(zdq).gt.prec) then
	534	zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq
	535	zsigg(ij,l)=1.-zsigd(ij,l)
	536	c if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
	537	c s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
	538	c print*,'probleme au point ij=',ij,' l=',l
	539	c print*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l)
	540	c print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
	541	c stop
	542	c endif
	543	else
	544	zsigd(ij,l)=0.5
	545	zsigg(ij,l)=0.5
	546	qd(ij,l)=q(ij,l)
	547	qg(ij,l)=q(ij,l)
	548	endif
	549	enddo
	550	enddo
	551
	552	c calcul de la pente maximum dans la maille en valeur absolue
	553
	554	do l=1,llm
	555	do ij=iip2,ip1jm-1
	556	if (u_m(ij,l).ge.0.) then
	557	zsigp=zsigd(ij,l)
	558	zsigm=zsigg(ij,l)
	559	zqp=qd(ij,l)
	560	zqm=qg(ij,l)
	561	zm=masse(ij,l)
	562	zq=q(ij,l)
	563	else
	564	zsigm=zsigd(ij+1,l)
	565	zsigp=zsigg(ij+1,l)
	566	zqm=qd(ij+1,l)
	567	zqp=qg(ij+1,l)
	568	zm=masse(ij+1,l)
	569	zq=q(ij+1,l)
	570	endif
	571	zu=abs(u_m(ij,l))
	572	ladvplus(ij,l)=zu.gt.zm
	573	zsig=zu/zm
	574	if(zsig.eq.0.) zsigp=0.1
	575	if (mode.eq.1) then
	576	if (zsig.le.zsigp) then
	577	u_mq(ij,l)=u_m(ij,l)*zqp
	578	else if (mode.eq.1) then
	579	u_mq(ij,l)=
	580	s sign(zm,u_m(ij,l))(zsigpzqp+(zsig-zsigp)*zqm)
	581	endif
	582	else
	583	if (zsig.le.zsigp) then
	584	u_mq(ij,l)=u_m(ij,l)(zqp-0.5zsig/zsigp*(zqp-zq))
	585	else
	586	zz=0.5*(zsig-zsigp)/zsigm
	587	u_mq(ij,l)=sign(zm,u_m(ij,l))( 0.5(zq+zqp)*zsigp
	588	s +(zsig-zsigp)(zq+zz(zqm-zq)) )
	589	endif
	590	endif
	591	c if(zsig.lt.0.) then
	592	c print*,'au point ij=',ij,' l=',l,' sig=',zsig
	593	c stop
	594	c endif
	595	enddo
	596	enddo
	597
	598	do l=1,llm
	599	do ij=iip1+iip1,ip1jm,iip1
	600	u_mq(ij,l)=u_mq(ij-iim,l)
	601	ladvplus(ij,l)=ladvplus(ij-iim,l)
	602	enddo
	603	enddo
	604
	605	c=================================================================
	606	C SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
	607	c=================================================================
	608	c tris des regions a traiter
	609	n0=0
	610	do l=1,llm
	611	nl(l)=0
	612	do ij=iip2,ip1jm
	613	if(ladvplus(ij,l)) then
	614	nl(l)=nl(l)+1
	615	u_mq(ij,l)=0.
	616	endif
	617	enddo
	618	n0=n0+nl(l)
	619	enddo
	620
	621	if(n0.gt.1) then
	622	IF (prt_level > 9) WRITE(lunout,*)
	623	& 'Nombre de points pour lesquels on advect plus que le'
	624	& ,'contenu de la maille : ',n0
	625
	626	do l=1,llm
	627	if(nl(l).gt.0) then
	628	iju=0
	629	c indicage des mailles concernees par le traitement special
	630	do ij=iip2,ip1jm
	631	if(ladvplus(ij,l).and.mod(ij,iip1).ne.0) then
	632	iju=iju+1
	633	indu(iju)=ij
	634	endif
	635	enddo
	636	niju=iju
	637	c print*,'niju,nl',niju,nl(l)
	638
	639	c traitement des mailles
	640	do iju=1,niju
	641	ij=indu(iju)
	642	j=(ij-1)/iip1+1
	643	zu_m=u_m(ij,l)
	644	u_mq(ij,l)=0.
	645	if(zu_m.gt.0.) then
	646	ijq=ij
	647	i=ijq-(j-1)*iip1
	648	c accumulation pour les mailles completements advectees
	649	do while(zu_m.gt.masse(ijq,l))
	650	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	651	zu_m=zu_m-masse(ijq,l)
	652	i=mod(i-2+iim,iim)+1
	653	ijq=(j-1)*iip1+i
	654	enddo
	655	c MODIFS SPECIFIQUES DU SCHEMA
	656	c ajout de la maille non completement advectee
	657	zsig=zu_m/masse(ijq,l)
	658	if(zsig.le.zsigd(ijq,l)) then
	659	u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l)
	660	s -0.5zsig/zsigd(ijq,l)(qd(ijq,l)-q(ijq,l)))
	661	else
	662	c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
	663	c goto 8888
	664	zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
	665	if(.not.(zz.gt.0..and.zz.le.0.5)) then
	666	WRITE(lunout,*)'probleme2 au point ij=',ij,
	667	s ' l=',l
	668	WRITE(lunout,*)'zz=',zz
	669	stop
	670	endif
	671	u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*(
	672	s 0.5(q(ijq,l)+qd(ijq,l))zsigd(ijq,l)
	673	s +(zsig-zsigd(ijq,l))(q(ijq,l)+zz(qg(ijq,l)-q(ijq,l))) )
	674	endif
	675	else
	676	ijq=ij+1
	677	i=ijq-(j-1)*iip1
	678	c accumulation pour les mailles completements advectees
	679	do while(-zu_m.gt.masse(ijq,l))
	680	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	681	zu_m=zu_m+masse(ijq,l)
	682	i=mod(i,iim)+1
	683	ijq=(j-1)*iip1+i
	684	enddo
	685	c ajout de la maille non completement advectee
	686	c 2eme MODIF SPECIFIQUE
	687	zsig=-zu_m/masse(ij+1,l)
	688	if(zsig.le.zsigg(ijq,l)) then
	689	u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l)
	690	s -0.5zsig/zsigg(ijq,l)(qg(ijq,l)-q(ijq,l)))
	691	else
	692	c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
	693	c goto 9999
	694	zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
	695	if(.not.(zz.gt.0..and.zz.le.0.5)) then
	696	WRITE(lunout,*)'probleme22 au point ij=',ij
	697	s ,' l=',l
	698	WRITE(lunout,*)'zz=',zz
	699	stop
	700	endif
	701	u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*(
	702	s 0.5(q(ijq,l)+qg(ijq,l))zsigg(ijq,l)
	703	s +(zsig-zsigg(ijq,l))*
	704	s (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
	705	endif
	706	c fin de la modif
	707	endif
	708	enddo
	709	endif
	710	enddo
	711	endif ! n0.gt.0
	712
	713	c bouclage en latitude
	714	do l=1,llm
	715	do ij=iip1+iip1,ip1jm,iip1
	716	u_mq(ij,l)=u_mq(ij-iim,l)
	717	enddo
	718	enddo
	719
	720	c=================================================================
	721	c CALCUL DE LA CONVERGENCE DES FLUX
	722	c=================================================================
	723
	724	do l=1,llm
	725	do ij=iip2+1,ip1jm
	726	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	727	q(ij,l)=(q(ij,l)*masse(ij,l)+
	728	& u_mq(ij-1,l)-u_mq(ij,l))
	729	& /new_m
	730	masse(ij,l)=new_m
	731	enddo
	732	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	733	do ij=iip1+iip1,ip1jm,iip1
	734	q(ij-iim,l)=q(ij,l)
	735	masse(ij-iim,l)=masse(ij,l)
	736	enddo
	737	enddo
	738
	739	RETURN
	740	END
	741	SUBROUTINE advny(q,qs,qn,masse,v_m)
	742	c
	743	c Auteur : F. Hourdin
	744	c
	745	c ********************************************************************
	746	c Shema d'advection " pseudo amont " .
	747	c ********************************************************************
	748	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	749	c
	750	c
	751	c --------------------------------------------------------------------
	752	IMPLICIT NONE
	753	c
	754	#include "dimensions.h"
	755	#include "paramet.h"
	756	#include "comgeom.h"
	757	#include "iniprint.h"
	758	c
	759	c
	760	c Arguments:
	761	c ----------
	762	real masse(ip1jmp1,llm)
	763	real v_m( ip1jm,llm )
	764	real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
	765	c
	766	c Local
	767	c ---------
	768	c
	769	INTEGER ij,l
	770	c
	771	real new_m,zdq,zz
	772	real zsigs(ip1jmp1),zsign(ip1jmp1),zsig
	773	real v_mq(ip1jm,llm)
	774	real convpn,convps,convmpn,convmps,massen,masses
	775	real zm,zq,zsigm,zsigp,zqm,zqp
	776	real ssum
	777	real prec
	778	save prec
	779
	780	#ifdef CRAY
	781	data prec/1.e-24/
	782	#else
	783	data prec/1.e-15/
	784	#endif
	785	do l=1,llm
	786	do ij=1,ip1jmp1
	787	zdq=qn(ij,l)-qs(ij,l)
	788	c if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
	789	c print*,'probleme au point ij=',ij,' l=',l,' advnqx'
	790	c print*,qn(ij,l),q(ij,l),qs(ij,l)
	791	c qn(ij,l)=q(ij,l)
	792	c qs(ij,l)=q(ij,l)
	793	c endif
	794	if(abs(zdq).gt.prec) then
	795	zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
	796	zsigs(ij)=1.-zsign(ij)
	797	c if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
	798	c s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
	799	c print*,'probleme au point ij=',ij,' l=',l
	800	c print*,'sigs=',zsigs(ij),' sign=',zsign(ij)
	801	c stop
	802	c endif
	803	else
	804	zsign(ij)=0.5
	805	zsigs(ij)=0.5
	806	endif
	807	enddo
	808
	809	c calcul de la pente maximum dans la maille en valeur absolue
	810
	811	do ij=1,ip1jm
	812	if (v_m(ij,l).ge.0.) then
	813	zsigp=zsign(ij+iip1)
	814	zsigm=zsigs(ij+iip1)
	815	zqp=qn(ij+iip1,l)
	816	zqm=qs(ij+iip1,l)
	817	zm=masse(ij+iip1,l)
	818	zq=q(ij+iip1,l)
	819	else
	820	zsigm=zsign(ij)
	821	zsigp=zsigs(ij)
	822	zqm=qn(ij,l)
	823	zqp=qs(ij,l)
	824	zm=masse(ij,l)
	825	zq=q(ij,l)
	826	endif
	827	zsig=abs(v_m(ij,l))/zm
	828	if(zsig.eq.0.) zsigp=0.1
	829	if (zsig.le.zsigp) then
	830	v_mq(ij,l)=v_m(ij,l)(zqp-0.5zsig/zsigp*(zqp-zq))
	831	else
	832	zz=0.5*(zsig-zsigp)/zsigm
	833	v_mq(ij,l)=sign(zm,v_m(ij,l))( 0.5(zq+zqp)*zsigp
	834	s +(zsig-zsigp)(zq+zz(zqm-zq)) )
	835	endif
	836	enddo
	837	enddo
	838
	839	do l=1,llm
	840	do ij=iip2,ip1jm
	841	new_m=masse(ij,l)
	842	& +v_m(ij,l)-v_m(ij-iip1,l)
	843	q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))
	844	& /new_m
	845	masse(ij,l)=new_m
	846	enddo
	847	c.-. ancienne version
	848	convpn=SSUM(iim,v_mq(1,l),1)
	849	convmpn=ssum(iim,v_m(1,l),1)
	850	massen=ssum(iim,masse(1,l),1)
	851	new_m=massen+convmpn
	852	q(1,l)=(q(1,l)*massen+convpn)/new_m
	853	do ij = 1,iip1
	854	q(ij,l)=q(1,l)
	855	masse(ij,l)=new_m*aire(ij)/apoln
	856	enddo
	857
	858	convps=-SSUM(iim,v_mq(ip1jm-iim,l),1)
	859	convmps=-ssum(iim,v_m(ip1jm-iim,l),1)
	860	masses=ssum(iim,masse(ip1jm+1,l),1)
	861	new_m=masses+convmps
	862	q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m
	863	do ij = ip1jm+1,ip1jmp1
	864	q(ij,l)=q(ip1jm+1,l)
	865	masse(ij,l)=new_m*aire(ij)/apols
	866	enddo
	867	enddo
	868
	869	RETURN
	870	END
	871	SUBROUTINE advnz(q,qh,qb,masse,w_m)
	872	c
	873	c Auteurs: F.Hourdin
	874	c
	875	c ********************************************************************
	876	c Shema d'advection " pseudo amont " .
	877	c b designe le bas et h le haut
	878	c il y a une correspondance entre le b en z et le d en x
	879	c ********************************************************************
	880	c
	881	c
	882	c --------------------------------------------------------------------
	883	IMPLICIT NONE
	884	c
	885	#include "dimensions.h"
	886	#include "paramet.h"
	887	#include "comgeom.h"
	888	#include "iniprint.h"
	889	c
	890	c
	891	c Arguments:
	892	c ----------
	893	real masse(ip1jmp1,llm)
	894	real w_m( ip1jmp1,llm+1)
	895	real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
	896
	897	c
	898	c Local
	899	c ---------
	900	c
	901	INTEGER ij,l
	902	c
	903	real new_m,zdq,zz
	904	real zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig
	905	real w_mq(ip1jmp1,llm+1)
	906	real zm,zq,zsigm,zsigp,zqm,zqp
	907	real prec
	908	save prec
	909
	910	#ifdef CRAY
	911	data prec/1.e-24/
	912	#else
	913	data prec/1.e-13/
	914	#endif
	915
	916	do l=1,llm
	917	do ij=1,ip1jmp1
	918	zdq=qb(ij,l)-qh(ij,l)
	919	c if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
	920	c print*,'probleme au point ij=',ij,' l=',l
	921	c print*,qh(ij,l),q(ij,l),qb(ij,l)
	922	c qh(ij,l)=q(ij,l)
	923	c qb(ij,l)=q(ij,l)
	924	c endif
	925
	926	if(abs(zdq).gt.prec) then
	927	zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
	928	zsigh(ij,l)=1.-zsigb(ij,l)
	929	zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)
	930	else
	931	zsigb(ij,l)=0.5
	932	zsigh(ij,l)=0.5
	933	endif
	934	enddo
	935	enddo
	936
	937	c print*,'ok1'
	938	c calcul de la pente maximum dans la maille en valeur absolue
	939	do l=2,llm
	940	do ij=1,ip1jmp1
	941	if (w_m(ij,l).ge.0.) then
	942	zsigp=zsigb(ij,l)
	943	zsigm=zsigh(ij,l)
	944	zqp=qb(ij,l)
	945	zqm=qh(ij,l)
	946	zm=masse(ij,l)
	947	zq=q(ij,l)
	948	else
	949	zsigm=zsigb(ij,l-1)
	950	zsigp=zsigh(ij,l-1)
	951	zqm=qb(ij,l-1)
	952	zqp=qh(ij,l-1)
	953	zm=masse(ij,l-1)
	954	zq=q(ij,l-1)
	955	endif
	956	zsig=abs(w_m(ij,l))/zm
	957	if(zsig.eq.0.) zsigp=0.1
	958	if (zsig.le.zsigp) then
	959	w_mq(ij,l)=w_m(ij,l)(zqp-0.5zsig/zsigp*(zqp-zq))
	960	else
	961	zz=0.5*(zsig-zsigp)/zsigm
	962	w_mq(ij,l)=sign(zm,w_m(ij,l))( 0.5(zq+zqp)*zsigp
	963	s +(zsig-zsigp)(zq+zz(zqm-zq)) )
	964	endif
	965	enddo
	966	enddo
	967
	968	do ij=1,ip1jmp1
	969	w_mq(ij,llm+1)=0.
	970	w_mq(ij,1)=0.
	971	enddo
	972
	973	do l=1,llm
	974	do ij=1,ip1jmp1
	975	new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)
	976	q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))
	977	& /new_m
	978	masse(ij,l)=new_m
	979	enddo
	980	enddo
	981	c print*,'ok3'
	982	RETURN
	983	END

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

Download in other formats:

Original Format