Context Navigation

advn.F @ 3189

Last change on this file since 3189 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.8 KB

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

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

Download in other formats:

Original Format