Context Navigation

vlspltqs.F @ 2421

Last change on this file since 2421 was 2286, checked in by Ehouarn Millour, 9 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
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 25.1 KB

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

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

Download in other formats:

Original Format