Context Navigation

vlspltqs.F @ 5

Last change on this file since 5 was 5, checked in by slebonnois, 14 years ago

Ajout des fichiers .def Venus et Titan (tels qu'ils sont utilisés

actuellement) dans les deftanks.

Ajout d'une doc sur Cp(T).
Modifications dans dyn3d concernant Cp(T), cf le log (v5) dans chantiers
Premières modifs de l'appel à la physique dans dyn3d/calfis, cf log (v5)
Elimination des cpdet.* dans phytitan et phyvenus (remplacés par cpdet.F

dans dyn3d).

File size: 21.1 KB

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

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

Download in other formats:

Original Format