Context Navigation

vlspltqs.F @ 3493

Last change on this file since 3493 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

File size: 25.0 KB

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

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

Download in other formats:

Original Format