Context Navigation

vlspltqs.F @ 4802

Last change on this file since 4802 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

Line
1	c
2	c $Id: vlspltqs.F 2286 2015-05-20 13:27:07Z lguez $
3	c
4	SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5	, p,pk,teta,iq )
6	USE infotrac, ONLY: nqtot,nqdesc,iqfils
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)
38	REAL q(ip1jmp1,llm,nqtot)
39	REAL w(ip1jmp1,llm),pdt
40	REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
41	INTEGER iq ! CRisi
42	c
43	c Local
44	c ---------
45	c
46	INTEGER i,ij,l,j,ii
47	INTEGER ifils,iq2 ! CRisi
48	c
49	REAL qsat(ip1jmp1,llm)
50	REAL zm(ip1jmp1,llm,nqtot)
51	REAL mu(ip1jmp1,llm)
52	REAL mv(ip1jm,llm)
53	REAL mw(ip1jmp1,llm+1)
54	REAL zq(ip1jmp1,llm,nqtot)
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
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
129
130	c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
131	call vlxqs(zq,pente_max,zm,mu,qsat,iq)
132
133	c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
134
135	call vlyqs(zq,pente_max,zm,mv,qsat,iq)
136
137	c call minmaxq(zq,qmin,qmax,'avant vlz ')
138
139	call vlz(zq,pente_max,zm,mw,iq)
140
141	c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
142	c call minmaxq(zm,qmin,qmax,'M avant vlyqs ')
143
144	call vlyqs(zq,pente_max,zm,mv,qsat,iq)
145
146	c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
147	c call minmaxq(zm,qmin,qmax,'M avant vlxqs ')
148
149	call vlxqs(zq,pente_max,zm,mu,qsat,iq)
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
157	q(ij,l,iq)=zq(ij,l,iq)
158	ENDDO
159	DO ij=1,ip1jm+1,iip1
160	q(ij+iim,l,iq)=q(ij,l,iq)
161	ENDDO
162	ENDDO
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
177	!write(,) 'vlspltqs 183: fin de la routine'
178
179	RETURN
180	END
181	SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
182	USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
183
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 ----------
203	REAL masse(ip1jmp1,llm,nqtot),pente_max
204	REAL u_m( ip1jmp1,llm )
205	REAL q(ip1jmp1,llm,nqtot)
206	REAL qsat(ip1jmp1,llm)
207	INTEGER iq ! CRisi
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
221	! CRisi
222	REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
223	INTEGER ifils,iq2 ! CRisi
224
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
256	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
257	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
258	c sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
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
310	dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
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
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),
354	, u_m(ij,l))
355	zdum(ij,l)=0.5*zdum(ij,l)
356	u_mq(ij,l)=cvmgp(
357	, q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
358	, q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
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
370	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
371	u_mq(ij,l)=u_m(ij,l)*
372	$ min(q(ij,l,iq)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
373	ELSE
374	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
375	u_mq(ij,l)=u_m(ij,l)*
376	$ min(q(ij+1,l,iq)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
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
418	IF(n0.gt.0) THEN
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
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)
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*
454	& (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
455	& *dxq(ijq,l))
456	ELSE
457	ijq=ij+1
458	i=ijq-(j-1)*iip1
459	c accumulation pour les mailles completements advectees
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)
464	i=mod(i,iim)+1
465	ijq=(j-1)*iip1+i
466	ENDDO
467	c ajout de la maille non completement advectee
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))
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
486	! CRisi: appel récursif de l'advection sur les fils.
487	! Il faut faire ça avant d'avoir mis à jour q et masse
488	!write(,) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
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
507
508	c calcul des tendances
509
510	DO l=1,llm
511	DO ij=iip2+1,ip1jm
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)+
514	& u_mq(ij-1,l)-u_mq(ij,l))
515	& /new_m
516	masse(ij,l,iq)=new_m
517	ENDDO
518	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
519	DO ij=iip1+iip1,ip1jm,iip1
520	q(ij-iim,l,iq)=q(ij,l,iq)
521	masse(ij-iim,l,iq)=masse(ij,l,iq)
522	ENDDO
523	ENDDO
524
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
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
548	SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
549	USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
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 ----------
573	REAL masse(ip1jmp1,llm,nqtot),pente_max
574	REAL masse_adv_v( ip1jm,llm)
575	REAL q(ip1jmp1,llm,nqtot)
576	REAL qsat(ip1jmp1,llm)
577	INTEGER iq ! CRisi
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
601
602	REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
603	INTEGER ifils,iq2 ! CRisi
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
642	airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
643	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
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
651	dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
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
666	dyq(ij,l)=qpns-q(ij+iip1,l,iq)
667	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
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)
709	C appn=abs(dyq(1)/dyqv(iip1+1))
710	C PRINT*,dyq(ip1jm+1)
711	C PRINT*,dyqv(ip1jm-iip1+1)
712	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
713	C DO ij=2,iim
714	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
715	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
716	C ENDDO
717	C appn=min(pente_max/appn,1.)
718	C apps=min(pente_max/apps,1.)
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.)
724	C & appn=0.
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.)
727	C & apps=0.
728	C
729	C limitation des pentes aux poles
730	C DO ij=1,iip1
731	C dyq(ij)=appn*dyq(ij)
732	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
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
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)))
788	ELSE
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)) )
791	ENDIF
792	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
793	ENDDO
794	ENDDO
795
796
797	! CRisi: appel récursif de l'advection sur les fils.
798	! Il faut faire ça avant d'avoir mis à jour q et masse
799	!write(,) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
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)
814	!write(,) 'vlyqs 783: appel rec de vly, iq2=',iq2
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
819	DO l=1,llm
820	DO ij=iip2,ip1jm
821	newmasse=masse(ij,l,iq)
822	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
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
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
831	newmasse=masse(ij,l,iq)+convmpn*aire(ij)
832	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpnaire(ij))/
833	& newmasse
834	masse(ij,l,iq)=newmasse
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
839	newmasse=masse(ij,l,iq)+convmps*aire(ij)
840	q(ij,l,iq)=(q(ij,l,iq)masse(ij,l,iq)+convpsaire(ij))/
841	& newmasse
842	masse(ij,l,iq)=newmasse
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
855	c masse(ij,l,iq)=newmasse*aire(ij)
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
865	c masse(ij,l,iq)=newmasse*aire(ij)
866	c ENDDO
867	c._. fin nouvelle version
868	ENDDO
869
870	!write(,) 'vly 866'
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
883	!write(,) 'vly 879'
884
885	RETURN
886	END

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

Download in other formats:

Original Format