Context Navigation

vlspltqs_loc.F @ 2396

Last change on this file since 2396 was 2298, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes -r2237:2291 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: 23.2 KB

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

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

Download in other formats:

Original Format