Context Navigation

vlspltqs.F @ 5456

Last change on this file since 5456 was 1520, checked in by Ehouarn Millour, 14 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 21.1 KB

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

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

Download in other formats:

Original Format