Context Navigation

advn.F90 @ 5273

Last change on this file since 5273 was 5272, checked in by abarral, 9 months ago
Turn paramet.h into a module
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: 24.8 KB

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

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

Download in other formats:

Original Format