Context Navigation

vlsplt_loc.F @ 2049

Last change on this file since 2049 was 1910, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r1860:1909 into testing branch
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: 24.2 KB

Rev	Line
[1707]	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 --------------------------------------------------------------------
[1864]	15	USE parallel_lmdz
[1632]	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 --------------------------------------------------------------------
[1864]	350	USE parallel_lmdz
[1632]	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
[1707]	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)
[1707]	548	C appn=abs(dyq(1)/dyqv(iip1+1))
[1632]	549	C PRINT*,dyq(ip1jm+1)
	550	C PRINT*,dyqv(ip1jm-iip1+1)
[1707]	551	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[1632]	552	C DO ij=2,iim
[1707]	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
[1707]	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.)
[1707]	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.)
[1707]	566	C & apps=0.
[1632]	567	C
	568	C limitation des pentes aux poles
	569	C DO ij=1,iip1
[1707]	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 --------------------------------------------------------------------
[1864]	739	USE parallel_lmdz
[1632]	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