Context Navigation

vlspltqs_loc.F @ 2353

Last change on this file since 2353 was 2286, checked in by Ehouarn Millour, 10 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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
[2270]	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 ********************************************************************
[2270]	6	c Shema d''advection " pseudo amont " .
[1632]	7	c ********************************************************************
	8	c
	9	c --------------------------------------------------------------------
[1823]	10	USE parallel_lmdz
[2270]	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 ----------
[2270]	23	REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
[1632]	24	REAL u_m( ijb_u:ije_u,llm )
[2270]	25	REAL q(ijb_u:ije_u,llm,nqtot)
[1632]	26	REAL qsat(ijb_u:ije_u,llm)
[2270]	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)
[2281]	40	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
[2270]	41	INTEGER ifils,iq2 ! CRisi
[1632]	42
[2270]	43
[1632]	44	REAL SSUM
	45
	46
	47	INTEGER ijb,ije,ijb_x,ije_x
	48
[2286]	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
[2270]	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'
[2270]	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
[2270]	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
[2270]	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(
[2270]	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.
[2270]	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
[2270]	209	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
[1632]	210	u_mq(ij,l)=u_m(ij,l)*
[2270]	211	$ min(q(ij,l,iq)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
[1632]	212	ELSE
[2270]	213	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
[1632]	214	u_mq(ij,l)=u_m(ij,l)*
[2270]	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
[2286]	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
[2270]	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
[2270]	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
[2270]	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
[2270]	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
[2270]	336	! CRisi: appel récursif de l'advection sur les fils.
	337	! Il faut faire ça avant d'avoir mis à jour q et masse
[2286]	338	!write(,) 'vlspltqs 336: iq,ijb_x,nqfils(iq)=',
	339	! & iq,ijb_x,nqfils(iq)
[2270]	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
[2281]	347	masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
[2270]	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)
[2286]	355	!write(,) 'vlxqs 349: on appelle vlx pour iq2=',iq2
[2281]	356	call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
[2270]	357	enddo !do ifils=1,nqfils(iq)
	358	endif !if (nqfils(iq).gt.0) then
	359	! end CRisi
	360
[2286]	361	!write(,) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x
[2270]	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
[2270]	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
[2270]	371	masse(ij,l,iq)=new_m
[1632]	372	ENDDO
[2270]	373	c Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous)
[1632]	374	DO ij=ijb+iip1-1,ije,iip1
[2270]	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
[2270]	380
[2286]	381	!write(,) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x
[2270]	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
[2286]	400	!write(,) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x
[2270]	401
[1632]	402	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
[2270]	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
[2270]	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 --------------------------------------------------------------------
[1823]	420	USE parallel_lmdz
[2270]	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 ----------
[2270]	434	REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
[1632]	435	REAL masse_adv_v( ijb_v:ije_v,llm)
[2270]	436	REAL q(ijb_u:ije_u,llm,nqtot)
[1632]	437	REAL qsat(ijb_u:ije_u,llm)
[2270]	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)
[2270]	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
[2281]	464	REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
[2270]	465	INTEGER ifils,iq2 ! CRisi
	466
[1632]	467	REAL SSUM
	468
	469	DATA first/.true./
	470	INTEGER ijb,ije
	471
[2270]	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
[2286]	479	!write(,) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
	480	! & ij,l,iq,ijb,q(ij,l,:)
[2270]	481	endif
	482
[1632]	483	IF(first) THEN
	484	PRINT*,'Shema Amont nouveau appele dans Vanleer '
[2270]	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
[2270]	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
[2270]	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
[2270]	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
[2270]	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
[2270]	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)
[1673]	626	C appn=abs(dyq(1)/dyqv(iip1+1))
[1632]	627	C PRINT*,dyq(ip1jm+1)
	628	C PRINT*,dyqv(ip1jm-iip1+1)
[1673]	629	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
[1632]	630	C DO ij=2,iim
[1673]	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
[1673]	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.)
[1673]	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.)
[1673]	644	C & apps=0.
[1632]	645	C
	646	C limitation des pentes aux poles
	647	C DO ij=1,iip1
[1673]	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
[2270]	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
[2270]	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
[2270]	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
[2270]	725	! CRisi: appel récursif de l'advection sur les fils.
	726	! Il faut faire ça avant d'avoir mis à jour q et masse
[2286]	727	!write(,) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
[2270]	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
[2281]	740	masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
[2270]	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)
[2281]	748	call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
[2270]	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
[2270]	763	newmasse=masse(ij,l,iq)
[1632]	764	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
[2270]	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
[2270]	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
[2270]	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
[2270]	779	masse(ij,l,iq)=newmasse
[1632]	780	ENDDO
	781
	782	ENDIF
	783
	784	IF (pole_sud) THEN
	785
[2270]	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
[2270]	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
[2270]	792	masse(ij,l,iq)=newmasse
[1632]	793	ENDDO
	794
	795	ENDIF
	796	c.-. fin ancienne version
	797
	798	c._. nouvelle version
[2270]	799	c convpn=SSUM(iim,qbyv(1,l,iq),1)
[1632]	800	c convmpn=ssum(iim,masse_adv_v(1,l),1)
[2270]	801	c oldmasse=ssum(iim,masse(1,l,iq),1)
[1632]	802	c newmasse=oldmasse+convmpn
[2270]	803	c newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
[1632]	804	c newmasse=newmasse/apoln
	805	c DO ij = 1,iip1
[2270]	806	c q(ij,l,iq)=newq
	807	c masse(ij,l,iq)=newmasse*aire(ij)
[1632]	808	c ENDDO
[2270]	809	c convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
[1632]	810	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
[2270]	811	c oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
[1632]	812	c newmasse=oldmasse+convmps
[2270]	813	c newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
[1632]	814	c newmasse=newmasse/apols
	815	c DO ij = ip1jm+1,ip1jmp1
[2270]	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
[2270]	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