Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

vlspltqs_loc.F

Last change on this file was 4469, checked in by Laurent Fairhead, 14 months ago
Replaced STOP instructions by calls to abort_gcm for a cleaner exit
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.3 KB

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

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

Download in other formats:

Original Format