Context Navigation

vlspltqs_loc.F @ 2228

Last change on this file since 2228 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

File size: 19.2 KB

Rev	Line
[1632]	1	SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
	2	c
	3	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	4	c
	5	c ********************************************************************
	6	c Shema d'advection " pseudo amont " .
	7	c ********************************************************************
	8	c
	9	c --------------------------------------------------------------------
[1823]	10	USE parallel_lmdz
[1632]	11	IMPLICIT NONE
	12	c
	13	#include "dimensions.h"
	14	#include "paramet.h"
	15	#include "logic.h"
	16	#include "comvert.h"
	17	#include "comconst.h"
	18	c
	19	c
	20	c Arguments:
	21	c ----------
	22	REAL masse(ijb_u:ije_u,llm),pente_max
	23	REAL u_m( ijb_u:ije_u,llm )
	24	REAL q(ijb_u:ije_u,llm)
	25	REAL qsat(ijb_u:ije_u,llm)
	26	c
	27	c Local
	28	c ---------
	29	c
	30	INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
	31	INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
	32	c
	33	REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
	34	REAL dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
	35	REAL zz(ijb_u:ije_u)
	36	REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
	37	REAL u_mq(ijb_u:ije_u,llm)
	38
	39	REAL SSUM
	40
	41
	42	INTEGER ijb,ije,ijb_x,ije_x
	43
	44
	45	c calcul de la pente a droite et a gauche de la maille
	46
	47	c ijb=ij_begin
	48	c ije=ij_end
	49
	50	ijb=ijb_x
	51	ije=ije_x
	52
	53	if (pole_nord.and.ijb==1) ijb=ijb+iip1
	54	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
	55
	56	IF (pente_max.gt.-1.e-5) THEN
	57	c IF (pente_max.gt.10) THEN
	58
	59	c calcul des pentes avec limitation, Van Leer scheme I:
	60	c -----------------------------------------------------
	61
	62	c calcul de la pente aux points u
	63
	64	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	65	DO l = 1, llm
	66	DO ij=ijb,ije-1
	67	dxqu(ij)=q(ij+1,l)-q(ij,l)
	68	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
	69	c sigu(ij)=u_m(ij,l)/masse(ij,l)
	70	ENDDO
	71	DO ij=ijb+iip1-1,ije,iip1
	72	dxqu(ij)=dxqu(ij-iim)
	73	c sigu(ij)=sigu(ij-iim)
	74	ENDDO
	75
	76	DO ij=ijb,ije
	77	adxqu(ij)=abs(dxqu(ij))
	78	ENDDO
	79
	80	c calcul de la pente maximum dans la maille en valeur absolue
	81
	82	DO ij=ijb+1,ije
	83	dxqmax(ij,l)=pente_max*
	84	, min(adxqu(ij-1),adxqu(ij))
	85	c limitation subtile
	86	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
	87
	88
	89	ENDDO
	90
	91	DO ij=ijb+iip1-1,ije,iip1
	92	dxqmax(ij-iim,l)=dxqmax(ij,l)
	93	ENDDO
	94
	95	DO ij=ijb+1,ije
	96	#ifdef CRAY
	97	dxq(ij,l)=
	98	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
	99	#else
	100	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
	101	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
	102	ELSE
	103	c extremum local
	104	dxq(ij,l)=0.
	105	ENDIF
	106	#endif
	107	dxq(ij,l)=0.5*dxq(ij,l)
	108	dxq(ij,l)=
	109	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
	110	ENDDO
	111
	112	ENDDO ! l=1,llm
	113	c$OMP END DO NOWAIT
	114
	115	ELSE ! (pente_max.lt.-1.e-5)
	116
	117	c Pentes produits:
	118	c ----------------
	119	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	120	DO l = 1, llm
	121	DO ij=ijb,ije-1
	122	dxqu(ij)=q(ij+1,l)-q(ij,l)
	123	ENDDO
	124	DO ij=ijb+iip1-1,ije,iip1
	125	dxqu(ij)=dxqu(ij-iim)
	126	ENDDO
	127
	128	DO ij=ijb+1,ije
	129	zz(ij)=dxqu(ij-1)*dxqu(ij)
	130	zz(ij)=zz(ij)+zz(ij)
	131	IF(zz(ij).gt.0) THEN
	132	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
	133	ELSE
	134	c extremum local
	135	dxq(ij,l)=0.
	136	ENDIF
	137	ENDDO
	138
	139	ENDDO
	140	c$OMP END DO NOWAIT
	141	ENDIF ! (pente_max.lt.-1.e-5)
	142
	143	c bouclage de la pente en iip1:
	144	c -----------------------------
	145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	146	DO l=1,llm
	147	DO ij=ijb+iip1-1,ije,iip1
	148	dxq(ij-iim,l)=dxq(ij,l)
	149	ENDDO
	150
	151	DO ij=ijb,ije
	152	iadvplus(ij,l)=0
	153	ENDDO
	154
	155	ENDDO
	156	c$OMP END DO NOWAIT
	157
	158	if (pole_nord) THEN
	159	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	160	DO l=1,llm
	161	iadvplus(1:iip1,l)=0
	162	ENDDO
	163	c$OMP END DO NOWAIT
	164	endif
	165
	166	if (pole_sud) THEN
	167	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	168	DO l=1,llm
	169	iadvplus(ip1jm+1:ip1jmp1,l)=0
	170	ENDDO
	171	c$OMP END DO NOWAIT
	172	endif
	173
	174	c calcul des flux a gauche et a droite
	175
	176	#ifdef CRAY
	177	c--pas encore modification sur Qsat
	178	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	179	DO l=1,llm
	180	DO ij=ijb,ije-1
	181	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
	182	, 1.+u_m(ij,l)/masse(ij+1,l),
	183	, u_m(ij,l))
	184	zdum(ij,l)=0.5*zdum(ij,l)
	185	u_mq(ij,l)=cvmgp(
	186	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
	187	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
	188	, u_m(ij,l))
	189	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
	190	ENDDO
	191	ENDDO
	192	c$OMP END DO NOWAIT
	193
	194	#else
	195	c on cumule le flux correspondant a toutes les mailles dont la masse
	196	c au travers de la paroi pENDant le pas de temps.
	197	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
	198	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	199	DO l=1,llm
	200	DO ij=ijb,ije-1
	201	IF (u_m(ij,l).gt.0.) THEN
	202	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
	203	u_mq(ij,l)=u_m(ij,l)*
	204	$ min(q(ij,l)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
	205	ELSE
	206	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
	207	u_mq(ij,l)=u_m(ij,l)*
	208	$ min(q(ij+1,l)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
	209	ENDIF
	210	ENDDO
	211	ENDDO
	212	c$OMP END DO NOWAIT
	213	#endif
	214
	215
	216	c detection des points ou on advecte plus que la masse de la
	217	c maille
	218	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	219	DO l=1,llm
	220	DO ij=ijb,ije-1
	221	IF(zdum(ij,l).lt.0) THEN
	222	iadvplus(ij,l)=1
	223	u_mq(ij,l)=0.
	224	ENDIF
	225	ENDDO
	226	ENDDO
	227	c$OMP END DO NOWAIT
	228
	229	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	230	DO l=1,llm
	231	DO ij=ijb+iip1-1,ije,iip1
	232	iadvplus(ij,l)=iadvplus(ij-iim,l)
	233	ENDDO
	234	ENDDO
	235	c$OMP END DO NOWAIT
	236
	237
	238
	239	c traitement special pour le cas ou on advecte en longitude plus que le
	240	c contenu de la maille.
	241	c cette partie est mal vectorisee.
	242
	243	c pas d'influence de la pression saturante (pour l'instant)
	244
	245	c calcul du nombre de maille sur lequel on advecte plus que la maille.
	246
	247	n0=0
	248	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	249	DO l=1,llm
	250	nl(l)=0
	251	DO ij=ijb,ije
	252	nl(l)=nl(l)+iadvplus(ij,l)
	253	ENDDO
	254	n0=n0+nl(l)
	255	ENDDO
	256	c$OMP END DO NOWAIT
	257
	258	cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
	259	cym IF(n0.gt.1) THEN
	260	cym IF(n0.gt.0) THEN
	261	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
	262	ccc & ,'contenu de la maille : ',n0
	263	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	264	DO l=1,llm
	265	IF(nl(l).gt.0) THEN
	266	iju=0
	267	c indicage des mailles concernees par le traitement special
	268	DO ij=ijb,ije
	269	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
	270	iju=iju+1
	271	indu(iju)=ij
	272	ENDIF
	273	ENDDO
	274	niju=iju
	275	c PRINT*,'niju,nl',niju,nl(l)
	276
	277	c traitement des mailles
	278	DO iju=1,niju
	279	ij=indu(iju)
	280	j=(ij-1)/iip1+1
	281	zu_m=u_m(ij,l)
	282	u_mq(ij,l)=0.
	283	IF(zu_m.gt.0.) THEN
	284	ijq=ij
	285	i=ijq-(j-1)*iip1
	286	c accumulation pour les mailles completements advectees
	287	do while(zu_m.gt.masse(ijq,l))
	288	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
	289	zu_m=zu_m-masse(ijq,l)
	290	i=mod(i-2+iim,iim)+1
	291	ijq=(j-1)*iip1+i
	292	ENDDO
	293	c ajout de la maille non completement advectee
	294	u_mq(ij,l)=u_mq(ij,l)+zu_m*
	295	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
	296	ELSE
	297	ijq=ij+1
	298	i=ijq-(j-1)*iip1
	299	c accumulation pour les mailles completements advectees
	300	do while(-zu_m.gt.masse(ijq,l))
	301	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
	302	zu_m=zu_m+masse(ijq,l)
	303	i=mod(i,iim)+1
	304	ijq=(j-1)*iip1+i
	305	ENDDO
	306	c ajout de la maille non completement advectee
	307	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
	308	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
	309	ENDIF
	310	ENDDO
	311	ENDIF
	312	ENDDO
	313	c$OMP END DO NOWAIT
	314	cym ENDIF ! n0.gt.0
	315
	316
	317
	318	c bouclage en latitude
	319	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	320	DO l=1,llm
	321	DO ij=ijb+iip1-1,ije,iip1
	322	u_mq(ij,l)=u_mq(ij-iim,l)
	323	ENDDO
	324	ENDDO
	325	c$OMP END DO NOWAIT
	326
	327	c calcul des tendances
	328	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	329	DO l=1,llm
	330	DO ij=ijb+1,ije
	331	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
	332	q(ij,l)=(q(ij,l)*masse(ij,l)+
	333	& u_mq(ij-1,l)-u_mq(ij,l))
	334	& /new_m
	335	masse(ij,l)=new_m
	336	ENDDO
	337	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
	338	DO ij=ijb+iip1-1,ije,iip1
	339	q(ij-iim,l)=q(ij,l)
	340	masse(ij-iim,l)=masse(ij,l)
	341	ENDDO
	342	ENDDO
	343	c$OMP END DO NOWAIT
	344	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
	345	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
	346
	347
	348	RETURN
	349	END
	350	SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
	351	c
	352	c Auteurs: P.Le Van, F.Hourdin, F.Forget
	353	c
	354	c ********************************************************************
	355	c Shema d'advection " pseudo amont " .
	356	c ********************************************************************
	357	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
	358	c qsat est un argument de sortie pour le s-pg ....
	359	c
	360	c
	361	c --------------------------------------------------------------------
[1823]	362	USE parallel_lmdz
[1632]	363	IMPLICIT NONE
	364	c
	365	#include "dimensions.h"
	366	#include "paramet.h"
	367	#include "logic.h"
	368	#include "comvert.h"
	369	#include "comconst.h"
	370	#include "comgeom.h"
	371	c
	372	c
	373	c Arguments:
	374	c ----------
	375	REAL masse(ijb_u:ije_u,llm),pente_max
	376	REAL masse_adv_v( ijb_v:ije_v,llm)
	377	REAL q(ijb_u:ije_u,llm)
	378	REAL qsat(ijb_u:ije_u,llm)
	379	c
	380	c Local
	381	c ---------
	382	c
	383	INTEGER i,ij,l
	384	c
	385	REAL airej2,airejjm,airescb(iim),airesch(iim)
	386	REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
	387	REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
	388	REAL qbyv(ijb_v:ije_v,llm)
	389
	390	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
	391	c REAL newq,oldmasse
	392	Logical first
	393	SAVE first
	394	c$OMP THREADPRIVATE(first)
	395	REAL convpn,convps,convmpn,convmps
	396	REAL sinlon(iip1),sinlondlon(iip1)
	397	REAL coslon(iip1),coslondlon(iip1)
	398	SAVE sinlon,coslon,sinlondlon,coslondlon
	399	SAVE airej2,airejjm
	400	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
	401	c$OMP THREADPRIVATE(airej2,airejjm)
	402	c
	403	c
	404	REAL SSUM
	405
	406	DATA first/.true./
	407	INTEGER ijb,ije
	408
	409	IF(first) THEN
	410	PRINT*,'Shema Amont nouveau appele dans Vanleer '
	411	first=.false.
	412	do i=2,iip1
	413	coslon(i)=cos(rlonv(i))
	414	sinlon(i)=sin(rlonv(i))
	415	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
	416	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
	417	ENDDO
	418	coslon(1)=coslon(iip1)
	419	coslondlon(1)=coslondlon(iip1)
	420	sinlon(1)=sinlon(iip1)
	421	sinlondlon(1)=sinlondlon(iip1)
	422	airej2 = SSUM( iim, aire(iip2), 1 )
	423	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
	424	ENDIF
	425
	426	c
	427
	428	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	429	DO l = 1, llm
	430	c
	431	c --------------------------------
	432	c CALCUL EN LATITUDE
	433	c --------------------------------
	434
	435	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
	436	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
	437	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
	438
	439	if (pole_nord) then
	440	DO i = 1, iim
	441	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
	442	ENDDO
	443	qpns = SSUM( iim, airescb ,1 ) / airej2
	444	endif
	445
	446	if (pole_sud) then
	447	DO i = 1, iim
	448	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
	449	ENDDO
	450	qpsn = SSUM( iim, airesch ,1 ) / airejjm
	451	endif
	452
	453
	454	c calcul des pentes aux points v
	455
	456	ijb=ij_begin-2*iip1
	457	ije=ij_end+iip1
	458	if (pole_nord) ijb=ij_begin
	459	if (pole_sud) ije=ij_end-iip1
	460
	461	DO ij=ijb,ije
	462	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
	463	adyqv(ij)=abs(dyqv(ij))
	464	ENDDO
	465
	466
	467	c calcul des pentes aux points scalaires
	468
	469	ijb=ij_begin-iip1
	470	ije=ij_end+iip1
	471	if (pole_nord) ijb=ij_begin+iip1
	472	if (pole_sud) ije=ij_end-iip1
	473
	474	DO ij=ijb,ije
	475	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
	476	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
	477	dyqmax(ij)=pente_max*dyqmax(ij)
	478	ENDDO
	479
	480	IF (pole_nord) THEN
	481
	482	c calcul des pentes aux poles
	483	DO ij=1,iip1
	484	dyq(ij,l)=qpns-q(ij+iip1,l)
	485	ENDDO
	486
	487	c filtrage de la derivee
	488	dyn1=0.
	489	dyn2=0.
	490	DO ij=1,iim
	491	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
	492	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
	493	ENDDO
	494	DO ij=1,iip1
	495	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
	496	ENDDO
	497
	498	c calcul des pentes limites aux poles
	499	fn=1.
	500	DO ij=1,iim
	501	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
	502	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
	503	ENDIF
	504	ENDDO
	505
	506	DO ij=1,iip1
	507	dyq(ij,l)=fn*dyq(ij,l)
	508	ENDDO
	509
	510	ENDIF
	511
	512	IF (pole_sud) THEN
	513
	514	DO ij=1,iip1
	515	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
	516	ENDDO
	517
	518	dys1=0.
	519	dys2=0.
	520
	521	DO ij=1,iim
	522	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
	523	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
	524	ENDDO
	525
	526	DO ij=1,iip1
	527	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
	528	ENDDO
	529
	530	c calcul des pentes limites aux poles
	531	fs=1.
	532	DO ij=1,iim
	533	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
	534	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
	535	ENDIF
	536	ENDDO
	537
	538	DO ij=1,iip1
	539	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
	540	ENDDO
	541
	542	ENDIF
	543
	544
	545	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	546	C En memoire de dIFferents tests sur la
	547	C limitation des pentes aux poles.
	548	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	549	C PRINT*,dyq(1)
	550	C PRINT*,dyqv(iip1+1)
[1673]	551	C appn=abs(dyq(1)/dyqv(iip1+1))
[1632]	552	C PRINT*,dyq(ip1jm+1)
	553	C PRINT*,dyqv(ip1jm-iip1+1)
[1673]	554	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[1632]	555	C DO ij=2,iim
[1673]	556	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
	557	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
[1632]	558	C ENDDO
[1673]	559	C appn=min(pente_max/appn,1.)
	560	C apps=min(pente_max/apps,1.)
[1632]	561	C
	562	C
	563	C cas ou on a un extremum au pole
	564	C
	565	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
[1673]	566	C & appn=0.
[1632]	567	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	568	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
[1673]	569	C & apps=0.
[1632]	570	C
	571	C limitation des pentes aux poles
	572	C DO ij=1,iip1
[1673]	573	C dyq(ij)=appn*dyq(ij)
	574	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
[1632]	575	C ENDDO
	576	C
	577	C test
	578	C DO ij=1,iip1
	579	C dyq(iip1+ij)=0.
	580	C dyq(ip1jm+ij-iip1)=0.
	581	C ENDDO
	582	C DO ij=1,ip1jmp1
	583	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
	584	C ENDDO
	585	C
	586	C changement 10 07 96
	587	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
	588	C & THEN
	589	C DO ij=1,iip1
	590	C dyqmax(ij)=0.
	591	C ENDDO
	592	C ELSE
	593	C DO ij=1,iip1
	594	C dyqmax(ij)=pente_max*abs(dyqv(ij))
	595	C ENDDO
	596	C ENDIF
	597	C
	598	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
	599	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
	600	C &THEN
	601	C DO ij=ip1jm+1,ip1jmp1
	602	C dyqmax(ij)=0.
	603	C ENDDO
	604	C ELSE
	605	C DO ij=ip1jm+1,ip1jmp1
	606	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
	607	C ENDDO
	608	C ENDIF
	609	C fin changement 10 07 96
	610	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	611
	612	c calcul des pentes limitees
	613	ijb=ij_begin-iip1
	614	ije=ij_end+iip1
	615	if (pole_nord) ijb=ij_begin+iip1
	616	if (pole_sud) ije=ij_end-iip1
	617
	618	DO ij=ijb,ije
	619	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
	620	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
	621	ELSE
	622	dyq(ij,l)=0.
	623	ENDIF
	624	ENDDO
	625
	626	ENDDO
	627	c$OMP END DO NOWAIT
	628
	629	ijb=ij_begin-iip1
	630	ije=ij_end
	631	if (pole_nord) ijb=ij_begin
	632	if (pole_sud) ije=ij_end-iip1
	633
	634	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	635	DO l=1,llm
	636	DO ij=ijb,ije
	637	IF( masse_adv_v(ij,l).GT.0. ) THEN
	638	qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) +
	639	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
	640	ELSE
	641	qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
	642	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
	643	ENDIF
	644	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
	645	ENDDO
	646	ENDDO
	647	c$OMP END DO NOWAIT
	648
	649	ijb=ij_begin
	650	ije=ij_end
	651	if (pole_nord) ijb=ij_begin+iip1
	652	if (pole_sud) ije=ij_end-iip1
	653
	654	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	655	DO l=1,llm
	656	DO ij=ijb,ije
	657	newmasse=masse(ij,l)
	658	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
	659	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
	660	& /newmasse
	661	masse(ij,l)=newmasse
	662	ENDDO
	663	c.-. ancienne version
	664
	665	IF (pole_nord) THEN
	666
	667	convpn=SSUM(iim,qbyv(1,l),1)/apoln
	668	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
	669	DO ij = 1,iip1
	670	newmasse=masse(ij,l)+convmpn*aire(ij)
	671	q(ij,l)=(q(ij,l)masse(ij,l)+convpnaire(ij))/
	672	& newmasse
	673	masse(ij,l)=newmasse
	674	ENDDO
	675
	676	ENDIF
	677
	678	IF (pole_sud) THEN
	679
	680	convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
	681	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
	682	DO ij = ip1jm+1,ip1jmp1
	683	newmasse=masse(ij,l)+convmps*aire(ij)
	684	q(ij,l)=(q(ij,l)masse(ij,l)+convpsaire(ij))/
	685	& newmasse
	686	masse(ij,l)=newmasse
	687	ENDDO
	688
	689	ENDIF
	690	c.-. fin ancienne version
	691
	692	c._. nouvelle version
	693	c convpn=SSUM(iim,qbyv(1,l),1)
	694	c convmpn=ssum(iim,masse_adv_v(1,l),1)
	695	c oldmasse=ssum(iim,masse(1,l),1)
	696	c newmasse=oldmasse+convmpn
	697	c newq=(q(1,l)*oldmasse+convpn)/newmasse
	698	c newmasse=newmasse/apoln
	699	c DO ij = 1,iip1
	700	c q(ij,l)=newq
	701	c masse(ij,l)=newmasse*aire(ij)
	702	c ENDDO
	703	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
	704	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
	705	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
	706	c newmasse=oldmasse+convmps
	707	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
	708	c newmasse=newmasse/apols
	709	c DO ij = ip1jm+1,ip1jmp1
	710	c q(ij,l)=newq
	711	c masse(ij,l)=newmasse*aire(ij)
	712	c ENDDO
	713	c._. fin nouvelle version
	714	ENDDO
	715	c$OMP END DO NOWAIT
	716	RETURN
	717	END

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

Download in other formats:

Original Format