Context Navigation

vlspltqs.F @ 1017

Last change on this file since 1017 was 1017, checked in by emillour, 11 years ago

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

File size: 21.1 KB

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

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

Download in other formats:

Original Format