Context Navigation

vlspltqs_loc.F @ 3958

Last change on this file since 3958 was 3800, checked in by Laurent Fairhead, 4 years ago
Modifications nécessaires pour les isotopes CRisi
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:keywords set to `Id`
File size: 24.7 KB

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

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

Download in other formats:

Original Format