Context Navigation

vlspltqs_p.F @ 1664

Last change on this file since 1664 was 1520, checked in by Ehouarn Millour, 14 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 24.6 KB

Rev	Line
[1520]	1	c
	2	c $Id: vlspltqs_p.F 1520 2011-05-23 11:37:09Z fairhead $
	3	c
[630]	4	SUBROUTINE vlspltqs_p ( q,pente_max,masse,w,pbaru,pbarv,pdt,
	5	, p,pk,teta )
	6	c
	7	c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
	8	c
	9	c ********************************************************************
	10	c Shema d'advection " pseudo amont " .
	11	c + test sur humidite specifique: Q advecte< Qsat aval
	12	c (F. Codron, 10/99)
	13	c ********************************************************************
	14	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
	15	c
	16	c pente_max facteur de limitation des pentes: 2 en general
	17	c 0 pour un schema amont
	18	c pbaru,pbarv,w flux de masse en u ,v ,w
	19	c pdt pas de temps
	20	c
	21	c teta temperature potentielle, p pression aux interfaces,
	22	c pk exner au milieu des couches necessaire pour calculer Qsat
	23	c --------------------------------------------------------------------
	24	USE parallel
	25	USE mod_hallo
	26	USE VAMPIR
	27	IMPLICIT NONE
	28
	29	c
	30	#include "dimensions.h"
	31	#include "paramet.h"
	32	#include "logic.h"
	33	#include "comvert.h"
	34	#include "comconst.h"
	35
	36	c
	37	c Arguments:
	38	c ----------
	39	REAL masse(ip1jmp1,llm),pente_max
	40	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
	41	REAL q(ip1jmp1,llm)
	42	REAL w(ip1jmp1,llm),pdt
	43	REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
	44	c
	45	c Local
	46	c ---------
	47	c
	48	INTEGER i,ij,l,j,ii
	49	c
	50	REAL qsat(ip1jmp1,llm)
	51	REAL zm(ip1jmp1,llm)
	52	REAL mu(ip1jmp1,llm)
	53	REAL mv(ip1jm,llm)
	54	REAL mw(ip1jmp1,llm+1)
	55	REAL zq(ip1jmp1,llm)
	56	REAL temps1,temps2,temps3
	57	REAL zzpbar, zzw
	58	LOGICAL testcpu
	59	SAVE testcpu
	60	SAVE temps1,temps2,temps3
	61
	62	REAL qmin,qmax
	63	DATA qmin,qmax/0.,1.e33/
	64	DATA testcpu/.false./
	65	DATA temps1,temps2,temps3/0.,0.,0./
	66
	67	c--pour rapport de melange saturant--
	68
	69	REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
	70	REAL ptarg,pdelarg,foeew,zdelta
	71	REAL tempe(ip1jmp1)
	72	INTEGER ijb,ije
	73	type(request) :: MyRequest1
	74	type(request) :: MyRequest2
	75
	76	c fonction psat(T)
	77
	78	FOEEW ( PTARG,PDELARG ) = EXP (
	79	* (R3LES(1.-PDELARG)+R3IESPDELARG) * (PTARG-RTT)
	80	* / (PTARG-(R4LES(1.-PDELARG)+R4IESPDELARG)) )
	81
	82	r2es = 380.11733
	83	r3les = 17.269
	84	r3ies = 21.875
	85	r4les = 35.86
	86	r4ies = 7.66
	87	retv = 0.6077667
	88	rtt = 273.16
	89
	90	c-- Calcul de Qsat en chaque point
	91	c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
	92	c pour eviter une exponentielle.
	93
	94	call SetTag(MyRequest1,100)
	95	call SetTag(MyRequest2,101)
	96
	97	ijb=ij_begin-iip1
	98	ije=ij_end+iip1
	99	if (pole_nord) ijb=ij_begin
	100	if (pole_sud) ije=ij_end
	101
	102
	103	DO l = 1, llm
	104	DO ij = ijb, ije
	105	tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
	106	ENDDO
	107	DO ij = ijb, ije
	108	zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
	109	play = 0.5*(p(ij,l)+p(ij,l+1))
	110	qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
	111	qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
	112	ENDDO
	113	ENDDO
	114
	115	c PRINT*,'Debut vlsplt version debug sans vlyqs'
	116
	117	zzpbar = 0.5 * pdt
	118	zzw = pdt
	119
	120	ijb=ij_begin
	121	ije=ij_end
	122	if (pole_nord) ijb=ijb+iip1
	123	if (pole_sud) ije=ije-iip1
	124
	125
	126	DO l=1,llm
	127	DO ij = ijb,ije
	128	mu(ij,l)=pbaru(ij,l) * zzpbar
	129	ENDDO
	130	ENDDO
	131
	132	ijb=ij_begin-iip1
	133	ije=ij_end
	134	if (pole_nord) ijb=ij_begin
	135	if (pole_sud) ije=ij_end-iip1
	136
	137	DO l=1,llm
	138	DO ij=ijb,ije
	139	mv(ij,l)=pbarv(ij,l) * zzpbar
	140	ENDDO
	141	ENDDO
	142
	143	ijb=ij_begin
	144	ije=ij_end
	145
	146	DO l=1,llm
	147	DO ij=ijb,ije
	148	mw(ij,l)=w(ij,l) * zzw
	149	ENDDO
	150	ENDDO
	151
	152	DO ij=ijb,ije
	153	mw(ij,llm+1)=0.
	154	ENDDO
	155
	156	c CALL SCOPY(ijp1llm,q,1,zq,1)
	157	c CALL SCOPY(ijp1llm,masse,1,zm,1)
	158
	159	ijb=ij_begin
	160	ije=ij_end
	161	zq(ijb:ije,1:llm)=q(ijb:ije,1:llm)
	162	zm(ijb:ije,1:llm)=masse(ijb:ije,1:llm)
	163
	164
	165	call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_begin,ij_begin+2*iip1-1)
	166	call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_end-2*iip1+1,ij_end)
	167
	168	call VTb(VTHallo)
	169	call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
	170	call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest1)
	171	call SendRequest(MyRequest1)
	172	call VTe(VTHallo)
	173
	174	call vlxqs_p(zq,pente_max,zm,mu,qsat,
	175	. ij_begin+2iip1,ij_end-2iip1)
	176
	177	call VTb(VTHallo)
	178	call WaitRecvRequest(MyRequest1)
	179	call VTe(VTHallo)
	180
	181	call vlyqs_p(zq,pente_max,zm,mv,qsat)
	182
	183	call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
	184	call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
	185
	186	call VTb(VTHallo)
	187	call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
	188	call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest2)
	189	call SendRequest(MyRequest2)
	190	call VTe(VTHallo)
	191
	192	call vlz_p(zq,pente_max,zm,mw,ij_begin+2iip1,ij_end-2iip1)
	193
	194	call VTb(VTHallo)
	195	call WaitRecvRequest(MyRequest2)
	196	call VTe(VTHallo)
	197
	198	call vlyqs_p(zq,pente_max,zm,mv,qsat)
	199
	200
	201	call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_begin,ij_end)
	202
	203
	204	ijb=ij_begin
	205	ije=ij_end
	206
	207	DO l=1,llm
	208	DO ij=ijb,ije
	209	q(ij,l)=zq(ij,l)
	210	ENDDO
	211	ENDDO
	212
	213	DO l=1,llm
	214	DO ij=ijb,ije-iip1+1,iip1
	215	q(ij+iim,l)=q(ij,l)
	216	ENDDO
	217	ENDDO
	218
	219	call WaitSendRequest(MyRequest1)
	220	call WaitSendRequest(MyRequest2)
	221
	222	RETURN
	223	END
	224	SUBROUTINE vlxqs_p(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
	225	c
	226	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	227	c
	228	c ********************************************************************
	229	c Shema d'advection " pseudo amont " .
	230	c ********************************************************************
	231	c
	232	c --------------------------------------------------------------------
	233	USE parallel
	234	IMPLICIT NONE
	235	c
	236	#include "dimensions.h"
	237	#include "paramet.h"
	238	#include "logic.h"
	239	#include "comvert.h"
	240	#include "comconst.h"
	241	c
	242	c
	243	c Arguments:
	244	c ----------
	245	REAL masse(ip1jmp1,llm),pente_max
	246	REAL u_m( ip1jmp1,llm )
	247	REAL q(ip1jmp1,llm)
	248	REAL qsat(ip1jmp1,llm)
	249	c
	250	c Local
	251	c ---------
	252	c
	253	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
	254	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
	255	c
	256	REAL new_m,zu_m,zdum(ip1jmp1,llm)
	257	REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
	258	REAL zz(ip1jmp1)
	259	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
	260	REAL u_mq(ip1jmp1,llm)
	261
	262	REAL SSUM
	263
	264
	265	INTEGER ijb,ije,ijb_x,ije_x
	266
	267
	268	c calcul de la pente a droite et a gauche de la maille
	269
	270	c ijb=ij_begin
	271	c ije=ij_end
	272
	273	ijb=ijb_x
	274	ije=ije_x
	275
	276	if (pole_nord.and.ijb==1) ijb=ijb+iip1
	277	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
	278
	279	IF (pente_max.gt.-1.e-5) THEN
	280	c IF (pente_max.gt.10) THEN
	281
	282	c calcul des pentes avec limitation, Van Leer scheme I:
	283	c -----------------------------------------------------
	284
	285	c calcul de la pente aux points u
[764]	286
	287	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	288	DO l = 1, llm
	289	DO ij=ijb,ije-1
	290	dxqu(ij)=q(ij+1,l)-q(ij,l)
	291	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	292	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	293	ENDDO
	294	DO ij=ijb+iip1-1,ije,iip1
	295	dxqu(ij)=dxqu(ij-iim)
	296	c sigu(ij)=sigu(ij-iim)
	297	ENDDO
	298
	299	DO ij=ijb,ije
	300	adxqu(ij)=abs(dxqu(ij))
	301	ENDDO
	302
	303	c calcul de la pente maximum dans la maille en valeur absolue
	304
	305	DO ij=ijb+1,ije
	306	dxqmax(ij,l)=pente_max*
	307	, min(adxqu(ij-1),adxqu(ij))
	308	c limitation subtile
	309	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	310
	311
	312	ENDDO
	313
	314	DO ij=ijb+iip1-1,ije,iip1
	315	dxqmax(ij-iim,l)=dxqmax(ij,l)
	316	ENDDO
	317
	318	DO ij=ijb+1,ije
	319	#ifdef CRAY
	320	dxq(ij,l)=
	321	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	322	#else
	323	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	324	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	325	ELSE
	326	c extremum local
	327	dxq(ij,l)=0.
	328	ENDIF
	329	#endif
	330	dxq(ij,l)=0.5*dxq(ij,l)
	331	dxq(ij,l)=
	332	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	333	ENDDO
	334
	335	ENDDO ! l=1,llm
[764]	336	c$OMP END DO NOWAIT
[630]	337
	338	ELSE ! (pente_max.lt.-1.e-5)
	339
	340	c Pentes produits:
	341	c ----------------
[764]	342	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	343	DO l = 1, llm
	344	DO ij=ijb,ije-1
	345	dxqu(ij)=q(ij+1,l)-q(ij,l)
	346	ENDDO
	347	DO ij=ijb+iip1-1,ije,iip1
	348	dxqu(ij)=dxqu(ij-iim)
	349	ENDDO
	350
	351	DO ij=ijb+1,ije
	352	zz(ij)=dxqu(ij-1)*dxqu(ij)
	353	zz(ij)=zz(ij)+zz(ij)
	354	IF(zz(ij).gt.0) THEN
	355	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	356	ELSE
	357	c extremum local
	358	dxq(ij,l)=0.
	359	ENDIF
	360	ENDDO
	361
	362	ENDDO
[764]	363	c$OMP END DO NOWAIT
[630]	364	ENDIF ! (pente_max.lt.-1.e-5)
	365
	366	c bouclage de la pente en iip1:
	367	c -----------------------------
[764]	368	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	369	DO l=1,llm
	370	DO ij=ijb+iip1-1,ije,iip1
	371	dxq(ij-iim,l)=dxq(ij,l)
	372	ENDDO
	373
	374	DO ij=ijb,ije
	375	iadvplus(ij,l)=0
	376	ENDDO
	377
	378	ENDDO
[764]	379	c$OMP END DO NOWAIT
[630]	380
[764]	381	if (pole_nord) THEN
	382	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	383	DO l=1,llm
	384	iadvplus(1:iip1,l)=0
	385	ENDDO
	386	c$OMP END DO NOWAIT
	387	endif
	388
	389	if (pole_sud) THEN
	390	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	391	DO l=1,llm
	392	iadvplus(ip1jm+1:ip1jmp1,l)=0
	393	ENDDO
	394	c$OMP END DO NOWAIT
	395	endif
	396
[630]	397	c calcul des flux a gauche et a droite
	398
	399	#ifdef CRAY
	400	c--pas encore modification sur Qsat
[764]	401	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	402	DO l=1,llm
	403	DO ij=ijb,ije-1
	404	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	405	, 1.+u_m(ij,l)/masse(ij+1,l),
	406	, u_m(ij,l))
	407	zdum(ij,l)=0.5*zdum(ij,l)
	408	u_mq(ij,l)=cvmgp(
	409	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	410	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	411	, u_m(ij,l))
	412	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	413	ENDDO
	414	ENDDO
[764]	415	c$OMP END DO NOWAIT
	416
[630]	417	#else
	418	c on cumule le flux correspondant a toutes les mailles dont la masse
	419	c au travers de la paroi pENDant le pas de temps.
	420	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
[764]	421	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	422	DO l=1,llm
	423	DO ij=ijb,ije-1
	424	IF (u_m(ij,l).gt.0.) THEN
	425	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	426	u_mq(ij,l)=u_m(ij,l)*
	427	$ min(q(ij,l)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
	428	ELSE
	429	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	430	u_mq(ij,l)=u_m(ij,l)*
	431	$ min(q(ij+1,l)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
	432	ENDIF
	433	ENDDO
	434	ENDDO
[764]	435	c$OMP END DO NOWAIT
[630]	436	#endif
	437
	438
	439	c detection des points ou on advecte plus que la masse de la
	440	c maille
[764]	441	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	442	DO l=1,llm
	443	DO ij=ijb,ije-1
	444	IF(zdum(ij,l).lt.0) THEN
	445	iadvplus(ij,l)=1
	446	u_mq(ij,l)=0.
	447	ENDIF
	448	ENDDO
	449	ENDDO
[764]	450	c$OMP END DO NOWAIT
	451
	452	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	453	DO l=1,llm
	454	DO ij=ijb+iip1-1,ije,iip1
	455	iadvplus(ij,l)=iadvplus(ij-iim,l)
	456	ENDDO
	457	ENDDO
[764]	458	c$OMP END DO NOWAIT
[630]	459
	460
	461
	462	c traitement special pour le cas ou on advecte en longitude plus que le
	463	c contenu de la maille.
	464	c cette partie est mal vectorisee.
	465
	466	c pas d'influence de la pression saturante (pour l'instant)
	467
	468	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	469
	470	n0=0
[764]	471	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	472	DO l=1,llm
	473	nl(l)=0
	474	DO ij=ijb,ije
	475	nl(l)=nl(l)+iadvplus(ij,l)
	476	ENDDO
	477	n0=n0+nl(l)
	478	ENDDO
[764]	479	c$OMP END DO NOWAIT
[630]	480
[764]	481	cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
[630]	482	cym IF(n0.gt.1) THEN
[764]	483	cym IF(n0.gt.0) THEN
[630]	484	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
	485	ccc & ,'contenu de la maille : ',n0
[764]	486	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	487	DO l=1,llm
	488	IF(nl(l).gt.0) THEN
	489	iju=0
	490	c indicage des mailles concernees par le traitement special
	491	DO ij=ijb,ije
	492	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	493	iju=iju+1
	494	indu(iju)=ij
	495	ENDIF
	496	ENDDO
	497	niju=iju
	498	c PRINT*,'niju,nl',niju,nl(l)
	499
	500	c traitement des mailles
	501	DO iju=1,niju
	502	ij=indu(iju)
	503	j=(ij-1)/iip1+1
	504	zu_m=u_m(ij,l)
	505	u_mq(ij,l)=0.
	506	IF(zu_m.gt.0.) THEN
	507	ijq=ij
	508	i=ijq-(j-1)*iip1
	509	c accumulation pour les mailles completements advectees
	510	do while(zu_m.gt.masse(ijq,l))
	511	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	512	zu_m=zu_m-masse(ijq,l)
	513	i=mod(i-2+iim,iim)+1
	514	ijq=(j-1)*iip1+i
	515	ENDDO
	516	c ajout de la maille non completement advectee
	517	u_mq(ij,l)=u_mq(ij,l)+zu_m*
	518	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
	519	ELSE
	520	ijq=ij+1
	521	i=ijq-(j-1)*iip1
	522	c accumulation pour les mailles completements advectees
	523	do while(-zu_m.gt.masse(ijq,l))
	524	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	525	zu_m=zu_m+masse(ijq,l)
	526	i=mod(i,iim)+1
	527	ijq=(j-1)*iip1+i
	528	ENDDO
	529	c ajout de la maille non completement advectee
	530	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
	531	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
	532	ENDIF
	533	ENDDO
	534	ENDIF
	535	ENDDO
[764]	536	c$OMP END DO NOWAIT
	537	cym ENDIF ! n0.gt.0
[630]	538
	539
	540
	541	c bouclage en latitude
[764]	542	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	543	DO l=1,llm
	544	DO ij=ijb+iip1-1,ije,iip1
	545	u_mq(ij,l)=u_mq(ij-iim,l)
	546	ENDDO
	547	ENDDO
[764]	548	c$OMP END DO NOWAIT
[630]	549
	550	c calcul des tendances
[764]	551	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	552	DO l=1,llm
	553	DO ij=ijb+1,ije
	554	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	555	q(ij,l)=(q(ij,l)*masse(ij,l)+
	556	& u_mq(ij-1,l)-u_mq(ij,l))
	557	& /new_m
	558	masse(ij,l)=new_m
	559	ENDDO
	560	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	561	DO ij=ijb+iip1-1,ije,iip1
	562	q(ij-iim,l)=q(ij,l)
	563	masse(ij-iim,l)=masse(ij,l)
	564	ENDDO
	565	ENDDO
[764]	566	c$OMP END DO NOWAIT
[630]	567	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	568	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	569
	570
	571	RETURN
	572	END
	573	SUBROUTINE vlyqs_p(q,pente_max,masse,masse_adv_v,qsat)
	574	c
	575	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	576	c
	577	c ********************************************************************
	578	c Shema d'advection " pseudo amont " .
	579	c ********************************************************************
	580	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	581	c qsat est un argument de sortie pour le s-pg ....
	582	c
	583	c
	584	c --------------------------------------------------------------------
	585	USE parallel
	586	IMPLICIT NONE
	587	c
	588	#include "dimensions.h"
	589	#include "paramet.h"
	590	#include "logic.h"
	591	#include "comvert.h"
	592	#include "comconst.h"
	593	#include "comgeom.h"
	594	c
	595	c
	596	c Arguments:
	597	c ----------
	598	REAL masse(ip1jmp1,llm),pente_max
	599	REAL masse_adv_v( ip1jm,llm)
	600	REAL q(ip1jmp1,llm)
	601	REAL qsat(ip1jmp1,llm)
	602	c
	603	c Local
	604	c ---------
	605	c
	606	INTEGER i,ij,l
	607	c
	608	REAL airej2,airejjm,airescb(iim),airesch(iim)
	609	REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
	610	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
	611	REAL qbyv(ip1jm,llm)
	612
	613	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	614	c REAL newq,oldmasse
[764]	615	Logical first
	616	SAVE first
	617	c$OMP THREADPRIVATE(first)
[630]	618	REAL convpn,convps,convmpn,convmps
	619	REAL sinlon(iip1),sinlondlon(iip1)
	620	REAL coslon(iip1),coslondlon(iip1)
	621	SAVE sinlon,coslon,sinlondlon,coslondlon
	622	SAVE airej2,airejjm
[764]	623	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
	624	c$OMP THREADPRIVATE(airej2,airejjm)
[630]	625	c
	626	c
	627	REAL SSUM
	628
[764]	629	DATA first/.true./
[630]	630	INTEGER ijb,ije
	631
	632	IF(first) THEN
	633	PRINT*,'Shema Amont nouveau appele dans Vanleer '
	634	first=.false.
	635	do i=2,iip1
	636	coslon(i)=cos(rlonv(i))
	637	sinlon(i)=sin(rlonv(i))
	638	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	639	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	640	ENDDO
	641	coslon(1)=coslon(iip1)
	642	coslondlon(1)=coslondlon(iip1)
	643	sinlon(1)=sinlon(iip1)
	644	sinlondlon(1)=sinlondlon(iip1)
	645	airej2 = SSUM( iim, aire(iip2), 1 )
	646	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
	647	ENDIF
	648
	649	c
	650
[764]	651	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	652	DO l = 1, llm
	653	c
	654	c --------------------------------
	655	c CALCUL EN LATITUDE
	656	c --------------------------------
	657
	658	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
	659	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
	660	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
	661
	662	if (pole_nord) then
	663	DO i = 1, iim
	664	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
	665	ENDDO
	666	qpns = SSUM( iim, airescb ,1 ) / airej2
	667	endif
	668
	669	if (pole_sud) then
	670	DO i = 1, iim
	671	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
	672	ENDDO
	673	qpsn = SSUM( iim, airesch ,1 ) / airejjm
	674	endif
	675
	676
	677	c calcul des pentes aux points v
	678
	679	ijb=ij_begin-2*iip1
	680	ije=ij_end+iip1
	681	if (pole_nord) ijb=ij_begin
	682	if (pole_sud) ije=ij_end-iip1
	683
	684	DO ij=ijb,ije
	685	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
	686	adyqv(ij)=abs(dyqv(ij))
	687	ENDDO
	688
	689
	690	c calcul des pentes aux points scalaires
	691
	692	ijb=ij_begin-iip1
	693	ije=ij_end+iip1
	694	if (pole_nord) ijb=ij_begin+iip1
	695	if (pole_sud) ije=ij_end-iip1
	696
	697	DO ij=ijb,ije
	698	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
	699	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
	700	dyqmax(ij)=pente_max*dyqmax(ij)
	701	ENDDO
	702
	703	IF (pole_nord) THEN
	704
	705	c calcul des pentes aux poles
	706	DO ij=1,iip1
	707	dyq(ij,l)=qpns-q(ij+iip1,l)
	708	ENDDO
	709
	710	c filtrage de la derivee
	711	dyn1=0.
	712	dyn2=0.
	713	DO ij=1,iim
	714	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
	715	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
	716	ENDDO
	717	DO ij=1,iip1
	718	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
	719	ENDDO
	720
	721	c calcul des pentes limites aux poles
	722	fn=1.
	723	DO ij=1,iim
	724	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
	725	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
	726	ENDIF
	727	ENDDO
	728
	729	DO ij=1,iip1
	730	dyq(ij,l)=fn*dyq(ij,l)
	731	ENDDO
	732
	733	ENDIF
	734
	735	IF (pole_sud) THEN
	736
	737	DO ij=1,iip1
	738	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
	739	ENDDO
	740
	741	dys1=0.
	742	dys2=0.
	743
	744	DO ij=1,iim
	745	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
	746	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
	747	ENDDO
	748
	749	DO ij=1,iip1
	750	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
	751	ENDDO
	752
	753	c calcul des pentes limites aux poles
	754	fs=1.
	755	DO ij=1,iim
	756	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
	757	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
	758	ENDIF
	759	ENDDO
	760
	761	DO ij=1,iip1
	762	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
	763	ENDDO
	764
	765	ENDIF
	766
	767
	768	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	769	C En memoire de dIFferents tests sur la
	770	C limitation des pentes aux poles.
	771	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	772	C PRINT*,dyq(1)
	773	C PRINT*,dyqv(iip1+1)
[1520]	774	C appn=abs(dyq(1)/dyqv(iip1+1))
[630]	775	C PRINT*,dyq(ip1jm+1)
	776	C PRINT*,dyqv(ip1jm-iip1+1)
[1520]	777	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[630]	778	C DO ij=2,iim
[1520]	779	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
	780	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
[630]	781	C ENDDO
[1520]	782	C appn=min(pente_max/appn,1.)
	783	C apps=min(pente_max/apps,1.)
[630]	784	C
	785	C
	786	C cas ou on a un extremum au pole
	787	C
	788	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
[1520]	789	C & appn=0.
[630]	790	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	791	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
[1520]	792	C & apps=0.
[630]	793	C
	794	C limitation des pentes aux poles
	795	C DO ij=1,iip1
[1520]	796	C dyq(ij)=appn*dyq(ij)
	797	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
[630]	798	C ENDDO
	799	C
	800	C test
	801	C DO ij=1,iip1
	802	C dyq(iip1+ij)=0.
	803	C dyq(ip1jm+ij-iip1)=0.
	804	C ENDDO
	805	C DO ij=1,ip1jmp1
	806	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
	807	C ENDDO
	808	C
	809	C changement 10 07 96
	810	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	811	C & THEN
	812	C DO ij=1,iip1
	813	C dyqmax(ij)=0.
	814	C ENDDO
	815	C ELSE
	816	C DO ij=1,iip1
	817	C dyqmax(ij)=pente_max*abs(dyqv(ij))
	818	C ENDDO
	819	C ENDIF
	820	C
	821	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	822	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	823	C &THEN
	824	C DO ij=ip1jm+1,ip1jmp1
	825	C dyqmax(ij)=0.
	826	C ENDDO
	827	C ELSE
	828	C DO ij=ip1jm+1,ip1jmp1
	829	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
	830	C ENDDO
	831	C ENDIF
	832	C fin changement 10 07 96
	833	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	834
	835	c calcul des pentes limitees
	836	ijb=ij_begin-iip1
	837	ije=ij_end+iip1
	838	if (pole_nord) ijb=ij_begin+iip1
	839	if (pole_sud) ije=ij_end-iip1
	840
	841	DO ij=ijb,ije
	842	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
	843	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
	844	ELSE
	845	dyq(ij,l)=0.
	846	ENDIF
	847	ENDDO
	848
	849	ENDDO
[764]	850	c$OMP END DO NOWAIT
[630]	851
	852	ijb=ij_begin-iip1
	853	ije=ij_end
	854	if (pole_nord) ijb=ij_begin
	855	if (pole_sud) ije=ij_end-iip1
	856
[764]	857	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	858	DO l=1,llm
	859	DO ij=ijb,ije
	860	IF( masse_adv_v(ij,l).GT.0. ) THEN
	861	qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) +
	862	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
	863	ELSE
	864	qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
	865	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
	866	ENDIF
	867	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
	868	ENDDO
	869	ENDDO
[764]	870	c$OMP END DO NOWAIT
[630]	871
	872	ijb=ij_begin
	873	ije=ij_end
	874	if (pole_nord) ijb=ij_begin+iip1
	875	if (pole_sud) ije=ij_end-iip1
[764]	876
	877	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]	878	DO l=1,llm
	879	DO ij=ijb,ije
	880	newmasse=masse(ij,l)
	881	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
	882	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
	883	& /newmasse
	884	masse(ij,l)=newmasse
	885	ENDDO
	886	c.-. ancienne version
	887
	888	IF (pole_nord) THEN
	889
	890	convpn=SSUM(iim,qbyv(1,l),1)/apoln
	891	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
	892	DO ij = 1,iip1
	893	newmasse=masse(ij,l)+convmpn*aire(ij)
	894	q(ij,l)=(q(ij,l)masse(ij,l)+convpnaire(ij))/
	895	& newmasse
	896	masse(ij,l)=newmasse
	897	ENDDO
	898
	899	ENDIF
	900
	901	IF (pole_sud) THEN
	902
	903	convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
	904	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
	905	DO ij = ip1jm+1,ip1jmp1
	906	newmasse=masse(ij,l)+convmps*aire(ij)
	907	q(ij,l)=(q(ij,l)masse(ij,l)+convpsaire(ij))/
	908	& newmasse
	909	masse(ij,l)=newmasse
	910	ENDDO
	911
	912	ENDIF
	913	c.-. fin ancienne version
	914
	915	c._. nouvelle version
	916	c convpn=SSUM(iim,qbyv(1,l),1)
	917	c convmpn=ssum(iim,masse_adv_v(1,l),1)
	918	c oldmasse=ssum(iim,masse(1,l),1)
	919	c newmasse=oldmasse+convmpn
	920	c newq=(q(1,l)*oldmasse+convpn)/newmasse
	921	c newmasse=newmasse/apoln
	922	c DO ij = 1,iip1
	923	c q(ij,l)=newq
	924	c masse(ij,l)=newmasse*aire(ij)
	925	c ENDDO
	926	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	927	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	928	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
	929	c newmasse=oldmasse+convmps
	930	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
	931	c newmasse=newmasse/apols
	932	c DO ij = ip1jm+1,ip1jmp1
	933	c q(ij,l)=newq
	934	c masse(ij,l)=newmasse*aire(ij)
	935	c ENDDO
	936	c._. fin nouvelle version
	937	ENDDO
[764]	938	c$OMP END DO NOWAIT
[630]	939	RETURN
	940	END

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

Download in other formats:

Original Format