Context Navigation

advn.F90 @ 5355

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

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

Download in other formats:

Original Format