Context Navigation

vlsplt_loc.F @ 1721

Last change on this file since 1721 was 1673, checked in by Laurent Fairhead, 12 years ago

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer

End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

File size: 24.1 KB

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

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

Download in other formats:

Original Format