Context Navigation

vlspltqs_p.F @ 5452

Last change on this file since 5452 was 2603, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn logic.h into module logic_mod.F90 EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 24.5 KB

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

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

Download in other formats:

Original Format