Context Navigation

vlspltqs.F @ 3020

Last change on this file since 3020 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

File size: 25.0 KB

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

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

Download in other formats:

Original Format